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STREAM-TUBE ANALYSIS USED TO CALCULATE 
ELECTRON CONCENTRATION IN THE PRESENCE OF ACCELERATING 
AND EVAPORATING WATER DROPLETS 

By Charles J. Schexnayder, Jr., John S. Evans, 
and Dennis M. Bushnell 
Langley Research Center 

SUMMARY 

An analysis is presented for use in the determination of the electron concentration 
in a flowing plasma in the presence of accelerating and evaporating water droplets. For 
given input gas -flow properties and gas composition, the resulting computer program is 
used to calculate gas density, gas and droplet velocity, temperature of the gas, enthalpy, 
fraction of droplet evaporated, Mach number, and gas composition as a function of distance 
along a stream tube. Representative calculations for the injection of water into a stream 
tube flowing over a reentry vehicle traveling at a velocity of 5.4 km/sec (17 717 ft/sec) 
are given. 


INTRODUCTION 

When a vehicle enters the atmosphere at high velocity, blackout of radio signals to 
or from the vehicle will occur due to free electrons in the plasma surrounding the vehicle 
(refs. 1 and 2). Many investigations in the past have dealt with methods of reducing the 
electron concentration in a plasma by means of material injection (refs. 3 to 13). One of 
the more promising ways of reducing the electron concentration is by the recombination 
of electrons and ions after the addition of water droplets into a plasma, as discussed in 
references 3, 5, 6, 8, and 12. This method works well in cool overionized plasmas where 
the concentration of free electrons is larger than it would be if there were equilibrium at 
the local gas temperature. 

The processes by which water droplets are placed in a flow field have been passed 
over because no way to include them in the analysis has been found. However, some 
understanding of these processes is necessary for practical application of water injection 
to relieve radio blackout. For instance, at the altitudes where blackout occurs, useful 
penetration distances into a supersonic flow can be achieved only by injection in the form 
of a high-velocity liquid jet, which is broken up into droplets by aerodynamic forces 
and/or boiling. Spatial distribution of the resultant droplets and the distribution of droplet 


sizes in a given stream tube are known to influence the effectiveness of water injection as 
an electron suppressant but are not available. However, experimental studies (ref. 14) 
of penetration and spreading permit the determination of average mass flux of water down- 
stream of the injection site. Measured values of mean droplet size (ref. 15) are available 
as a substitute for knowledge of droplet size distribution. 

A computer program has been written by Barbara L. Weigel of the Langley Research 
Center to implement a method for calculating the properties of a one -dimensional stream- 
tube flow of a mixture of air and water droplets. (See the appendix.) The program takes 
into consideration acceleration and evaporation of the droplets along with a chemical 
kinetic reaction system based on 15 chemical species and 40 homogeneous chemical reac- 
tions. In addition, five heterogeneous reactions are included to account for the recom- 
bination of electrons and ions at droplet surfaces. The program is used to compute gas 
density, velocity, temperature, composition, percent evaporation of water, and the velocity 
of droplets as a function of distance from the injection site. A derivation of the necessary 
equations is included. A program listing with a table of the necessary inputs is also 
given. Some typical results are presented. 

SYMBOLS 


A 

A 

Bl 

b 2 

B3 

B ij 

C D 

C D,C 


stream-tube area, cm 2 


constant in drag-coefficient expression (eq. (51)) 


16 


mass fraction of oxygen, - 
mass fraction of nitrogen, 
mass fraction of hydrogen, 


(y 1 + 2 y 2 + r 5 + r 6 + r 8 + r 10 + 713), 


16 o * 
+ 18 fW 


1 +fw* 

14 ( 2y 3 + ^5 + r 7 + y 10 + 2 ^11 + ^12), 


1 + fw 


2 » * 

I8 tw 

1 + fw* 


term representing contribution of jth reaction to disappearance of ith species, 
mol 2 /g2 

droplet drag coefficient 
continuum droplet drag coefficient 


Cp free molecular droplet drag coefficient 
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( G d)m= 0 droplet drag coefficient for M = 0 
Cp heat capacity at constant pressure, ergs/mol-°K 

Cp specific heat at constant pressure, ergs/g-°K 

Cp specific heat of mixture at constant pressure, ergs/g-°K 

D reaction energy (divided by RTq), dimensionless 


3 N 


Fi 

f 

fr 

f r,0 


G ii 


reference length, given as nose diameter, cm 


Ej_, . . ., Eg correction factors applied at end of each integration step to preserve 
correct mass fraction of total oxygen, nitrogen, and hydrogen 

F force on droplet, dynes 

F e efficiency factor for capture of electrons by water droplet 


free energy of ith species, ergs/mol 
mass fraction of water evaporated, 1 - r 3 
free molecular recovery factor 
droplet recovery factor for zero mass transfer 

term representing contribution of jth reaction to production of ith species, 

molVg 2 


g plasma-sheath thickness transition function, g b + 

g a = 3.8738 + 0.03432U + 0.19754U 2 - 0.10445U 3 + 0.013263U 4 
g b = 5.46466 - 0.09100U - 0.33579U 2 + 0.11571U 3 - 0.012074U 4 


s a -Sb 


1 + (£) 




Hi 


enthalpy of ith species, ergs/mol 
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h 


specific enthalpy, ergs/g 


h* 

i 

i 

K 

k 

L 

M 

M 

(MW)j 

m 

n A 

N e 

N Nu 

N Nu,C 

N Nu,FM 

%e 

N St 

n 

P 
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total specific enthalpy of shocked gas, ergs/g 
species identification number 
reaction identification number 
equilibrium constant, mol/cm^ or dimensionless 

specific reaction rate constant, cm^/mol-sec; also, thermal conductivity, 
ergs/cm-sec-°K 

latent heat of water, ergs/g 

Mach number of gas relative to droplet; also, third body in chemical reactions 
mean molecular weight of ions, g/mol 
molecular weight of species i, g/mol 
mass-flow rate, g/sec 

Avogadro constant, 6.02486 x 10^ particles/mol 
electron concentration, electrons/cm^ 

Nusselt number for heat transfer to droplet 
continuum Nusselt number 
free molecular Nusselt number 

effective Reynolds number, £ ( u , t listed as RES in computer program 

M i 

free molecular Stanton number; listed as ST in computer program 
constant in drag -coefficient expression (eq. (51)) 
pressure, dynes/cm^ 



q 

R 

r 

r 

S 


s 

T 

t 

U 


u 

V 


v 


Wi 


X 

x 

Z 


“ij 

r i 


heat -transfer rate to droplet, ergs/cm^-sec 

3 

universal gas constant, 8.31696 X 10^ ergs/mol-°K or 82.1023 g tm-cm 

°K 

droplet radius, cm 
droplet radius divided by rg 
speed ratio, — li ~ ' 

j 2RT I n 

distance along stream tube from bow shock on vehicle, cm 
temperature, °K 
time, sec 

normalized relative velocity between gas and droplet 
gas velocity, cm/sec 
velocity, cm/sec 
droplet velocity, cm/sec 

net rate of production or disappearance of species i, mol/cm^-sec 
initial ratio of water mass flow to gas mass flow 
mole fraction 

S " S 2 

distance along streamline path from injection site divided by d^, — - 

d N 

compressibility factor 

factor which selects nonzero contributions to Wj 
specific-heat ratio of ith species, 


( c pA 


_( C P/ R )i 

- l" 


jl 
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y ^ specific concentration of ith species, mol/g 

pp 

X n Debye length, 6.901 -±-, cm 

u l|Me 

[i reduced mass of colliding species, g/mol; also, water viscosity, 

dynes-sec/cm^ 

p density, g/cm^ 

ct reaction rate cross section, cm^ 

Q, collision integral for viscosity 

Subscripts: 

av average 

aw zero heat transfer at droplet surface 

f evaluated at mean temperature of droplet film 

1 water 

m gas mixture after injection 

rel relative 

s saturated vapor 

v vapor 

w droplet surface 

oo free-stream condition 

2 value at injection site 

0 zero mass transfer; also, standard temperature 

Dot over symbol denotes derivative with respect to time. 

Arrow over symbol denotes vector. 
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METHOD OF CALCULATION 

The equations derived herein have been programed in FORTRAN IV language and 
a description of the program is given in the appendix. As the program is presently 
written, the following limitations should be considered when it is being used: 

(1) Although the liquid specified for injection is water, any suitable liquid could be 
used as long as the liquid droplet properties are known and the necessary changes are 
made in the finite-rate chemistry system. 

(2) The program cannot be used with injection of liquids at the stagnation point in 
the flow. 

(3) The finite-rate chemistry system does not take into account negative ions, such 
as H~, OH”, O", O 2 ", NO", and NC> 2 ~, but is otherwise believed to be valid up to speeds 
of about 9.14 km/sec (30 000 ft /sec). 

(4) The use of a stream -tube analysis limits application to inviscid flows. 

Thus, the growth of the boundary layer as altitude increases must be taken into account. 

Conservation Equations 

By passing over the difficult and complex questions which surround the actual injec- 
tion of water into the flow field, the problem is reduced to the properties along a stream 
tube in which, at some initial station, water droplets have been deposited by an atomizing 
liquid jet. Due to the complicated two-phase nature of the problem, it is necessary to 
make several idealizations, among which are 

(1) The flow is steady, one-dimensional, and inviscid. 

(2) The droplet size, velocity, and distribution at any streamwise station are uniform. 

(3) Nusselt number and drag coefficient for the droplets are the same as for an 
isolated sphere. 

(4) The ratio of mass flow of water (liquid and vapor) to mass flow of air in the 
stream tube is constant; that is, downstream of the initial station no droplets enter or 
leave the stream tube. 

(5) The gas phase is a uniform mixture of ideal gases. 

(6) The pressure distribution along the stream tube is assumed to be known (see, 
for example, fig. 1). 

Analysis of the stream-tube flow begins at the point where water droplets are sud- 
denly assumed to be present. If injection-induced perturbations to the properties of a 
flow field are to be included, they must appear in the influence they have on the initial 
conditions specified at the injection station. 
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A control volume approach is used in which the integral form of the conservation 
laws are applied to a one-dimensional stream-tube flow of a mixture of liquid droplets and 
gas. (See sketch a). (Information concerning derivations similar to those that follow is 
available in refs. 5, 16, and 17, and these references should be consulted for complete 



station -^2 A 

Sketch a 

Conservation of mass.- The conservation of mass for steady flow can be defined as 




• dA = 0 


or 


m l, 2 + m m,2 = m m + m i 

Inasmuch as droplets are assumed to be of a uniform size, 


-M = 1 - f = 


m 7 


m 


l l, 2 


By defining 


m 7 o 

w * = 

*m,2 

equation (1) becomes 

m T 


1 + w* = * m 1 /1 f '”’* 


+ (1 - f)w 


m 


m,2 


or 


m 


“_ = i +fw * 


m 


m,2 


( 1 ) 

( 2 ) 


( 3 ) 


( 4 ) 


and, therefore, 

A 2 pu 

where the volume of the droplets is negligible as compared with the volume of the air- 
water mixture. 

Conservation of momentum.- The conservation of momentum for steady flow can be 
expressed in the form 

F = j VpV • dA 


or 


m m,2 u 2 + m Z,2 v 2 + 


/Po + P\ 

p 2 A 2 + [— — j( A ' A 2) = PA +m m u + (1 - f)m^ 2 v 


(5) 


In equation (5) the option of specifying pressure or area is given. For the present work 
the pressure distribution is specified. Dividing equation (5) by m m 2 = P 2 U 2 A 2 an< 3 


using equation (4) yield 
P 2 -P 


u = 


2pu 


— P + u 2 + fv 2 - (1 - f)v]w*l 

1 + fw* \2p 2 u 2 2 L 2 J f 


( 6 ) 


where p = p(x) is assumed to be known. 

Conservation of energy.- The conservation of energy for steady flow can be defined 


as 


or 




m. 


W 2 


+ h pV • dA = 0 


+ h- 


I Ur 


i, 2) + “m,2^-f- + h 2 ) = ™m(jr + h ) + <* " + M™*, 2 ( 7 ) 


It has been assumed in equation (7) that the average temperature of the liquid droplet does 
not change appreciably; that is, h^ = h^ 2 . Also, it has been assumed that all the heat 
input to the droplet goes into evaporation. By solving equation (7) for h, using equa- 
tion (3), and assuming that v 2 is small enough to be neglected, one obtains 


where 


h = 


h* + hjfw* 




2 * 

vw* 2 


1 + fw* 


(8) 


h, = h 


- L 


(9) 
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Initially the droplet velocity is computed from the following relationship: 


v 2 = 


600 p 2 u 2 w* 

7r(r) 3 p 


+ 1 


(10) 


Chemical Relationships 

The properties of the gas mixture after water injection are given by the equations 
that follow. 


For a mixture of perfect gases, the equation of state is 
P 


P = 


RT 


15 


(ID 


i=l 

The enthalpy of a mixture of ideal gases is defined as 

15 

'Hi'' 


h = RT 


hi 

i=l 


RT 


(12) 


where H^/RT is obtained by fitting the data of references 18 and 19 to a set of fourth- 
order polynomial equations for enthalpy, entropy, heat capacity, and free energy as 
described in reference 20. 

The rate equations governing the change in chemical composition have the form 
d^i d N 

— - = — Wj (13) 

dx pu 1 


where 




W 




(14) 


and 


w « - “ijA( B U - G *i) 


(15) 


The factor is equal to 1 if reaction j involves species i is equal to 2 for 

atoms produced in dissociation of C> 2 , N 2 , and H 2 ); otherwise, Qf^ is equal to 0. The 
species are listed in table I. The specific reaction rate constants kj are obtained from 
data presented in table n. Reactions written in this table are given as forward rates in 
the endothermic direction. The constant kj for the first reaction given in table II can 
be expressed in the form (ref. 21). 
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I 



(16) 


k j = N A ct ] 


8RT 

m 


1/2 




exp 


D i T 0 


and for reactions 2 through 6 in the form (ref. 21) 

if 72 " 1 


/8RT\ 1//2 

4e -1 j 

w / 



+ V J 


exp - 


0 

T 


(17) 


On the basis of a simple collision model, the following form (ref. 22) holds for reactions 7 
through 40: 

\l/2 ( DjT 0 \ 


“ p 


(18) 


The ce values for the preceding reactions represent the best available values based on 
information in references 21 and 22 and were adjusted for experimental data if any were 
available. The Dj values were obtained from reference 23. In addition to information 
in references 21 and 22 some helpful summaries and discussions on reaction rates can be 
found in references 24 to 26. The set of reactions in table n is for an air-water mixture 
flowing over a nonablating body. 

For the dissociation and recombination of O2, the Bjj and Gjj terms in equa- 
tion (15) have the following form: 


py 8 Z y i 


-b 8 2 = K2 

(19) 

-g 8 2 = >'2 X 'n 

(20) 


2 2 _ 

For the formation and disappearance of electrons by reaction 16, they have the form 

y n y is 


B 7 16 = _B 15 16 


K 


16 


G 7 16 - -G 15 16 - y 7 y 7 


( 21 ) 

( 22 ) 


The equilibrium constants for the reactions in table II are computed from the free energy 
of the species as defined by the previously mentioned polynomials and have the form 
similar to the following example: 


K 2 = exp 


f_£L 

.El 

h £2\ 

( 1 /l 

\ RT 

RT 

RT/ 

\82.1023T/ 


(23) 
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In addition to the 40 homogeneous gas reactions, the recombination of ions with 
electrons on the surfaces of water droplets can be included in the reaction system by 
defining net rates in the form 

Wy = P 2 k jri (24) 

where j refers to the five catalytic reactions and i refers to the ions NO + , N 2 + , N + , 
0 + , and H + . The rate constants for reactions 41 to 45 can be expressed in the form 

k; = 3ur 2 w* L 22 x loS^TFo) (25) 

3 4r 2 p^v(l + fw*J ' ' 

where F e (refs. 3 and 6) can either be a constant or can be computed as follows: 

F e = e" g 

where 


g = g b 


+ 


g a " g b 



(26) 


The exponent g is proportional to the electric potential attained by a droplet in the 
plasma. This potential has different limiting values (represented by g & and g b ) 
depending upon the ratio of plasma-sheath thickness to drop radius (ref. 6), and the expres- 
sion for g provides a convenient (but arbitrary) transition between the two limits. The 
limits g a and g b are themselves functions of the velocity of the droplets relative to the 
gas. Therefore, g a and g b are fitted with fourth-order polynomials as follows: 

g a = 3.8738 + 0.03432U + 0.19754U 2 - 0.10445U 3 + 0.013263U 4 (0£U§3) (27) 

g b = 5.46466 - 0.091000U - 0.33579U 2 + 0.11571U 3 - 0.012074U 4 (0 2 u g 3) (28) 

where 


U = 




(29) 


D 



(30) 
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14 


M = 


r i (MW) i 


i=10 


14 


n 


i=10 


(31) 


In order to ensure that the correct mass fractions of total oxygen, nitrogen, and 
hydrogen were maintained in the mixture, y i values were replaced at the end of each 
computational step with 


E 


m = 


16 


ll 

(Vi + 2 y 2 + r 5 + yq + yg + r io + y i3) 


(i = 2, 8, 13) (32) 


E 2n 


Br 


!4(2y 3 + y 5 + y 7 + y 1Q + 2 y u + y 12 ) 




(i = 3, 7, 11, 12) (33) 


E 3n = 


2 y 1 + 2y 4 + y 6 + y 9 + r 14 




(i = 4, 9, 14) (34) 


E 4 y t = 



+ lGE]^ 


18 


(i = 1) (35) 


E 5n = 


14E2 + 16Ej 
30 


n 


(i = 5, 10) (36) 


E 6ri = 


E3 + 16Ej 
17 


n 


When these corrections are made, 
15 

^ yi(MW)j = 1 

i=l 

is satisfied. 


(i = 6) (37) 


The preceding equations apply for a mixture of air, water vapor, and the products of 
reactions involving both. They are used in the program to calculate the changes which 
occur in the gaseous mixture during a forward step. However, dilution of the mixture by 
evaporation during a step is accounted for only at the end of each step. At the end of a step 
a new value of f is computed from equation (2). This f value leads to new values of 
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Bj, B 2 , and B 3 and also to an increase in the concentration variable for water vapor 
7\ given by + 

( y l)new = ( 7l ^old + SK f new “ f old) ( 38 ) 

Because of the revised values of Bj, B 2 , B 3 , and it is necessary to recompute the 
correction factors which ensure conservation of total oxygen, nitrogen, and hydrogen. 
After these calculations are made, the program proceeds with computations for the next 
step. 


Conversion of specific concentration to mole fraction can be achieved by using 


X = — — 


15 


(39) 


u 


i=l 


and electron number density can be obtained from 

N e = N A py 15 (40) 

The compressibility factor for the air and water mixture can be computed by the fol- 
lowing relationship: 

15 

u 


z = 


i=l 


!(>■*)..♦ Ml* 

where + (y^j - 0.034674063. 


(l + fw*) 


(41) 


Auxiliary Relationships 


Droplet accelera tion.- The equation of motion for a mean droplet (based on ref. 5) 
assumed to be spherical in shape is (for an isolated spherical droplet) 

1 o o 4uT 3 p,V 


or 




dv = 3(u - v) pC D d N 
dx 8 p^vr 


(v* 0 ) 


with r = rgf. The drag coefficient (ref. 5) was taken for M < 0.5 as 

t r \ . 51. 1M 

( Cd )m=0 + 

C D = 


N- 


Re 


1 + 0.256M 


( c °)m=o 


51. 1M 


N- 


Re 


(42) 


(43) 


14 


III 


with 




In order to use equation (48) , one needs 



(obtained from ref. 27) and 



The specific heat of the mixture is 



For M § 0.5 the drag coefficient (ref. 5) is given by 

C D = C D,C + ( C D,FM " C D,c) exp (“ AN Re) 
where C D ^, A, and n are given as a function of M in table HE. 


(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 
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Heat transfer to droplets .- In accordance with the discussion of the global energy 
equation, it is assumed that all the heat transferred to the droplet goes into evaporation 
and, thus, 



dr 3 

dt 


or 

dr 2 _ 2 q rd N 
dx P/^v 

By defining an average droplet Nusselt number as 

2qr5 p 

%u = k f (h aw - h w ) 

the following equation can be used: 

dr 2 _ N Nukf d N( h aw ~ h w) 
dx Lp^CpV 


(52) 


(53) 


(54) 


In the present case, 


Vi h - f ( u ~ v ) + c T - c T (55) 

h aw “ h w ~ t r,0 ^-4-Cpi c p i w 

which is the enthalpy available for droplet evaporation. By assuming that the Prandtl 
number is unity, the following relation can be utilized: 


itf - h f c pjf 


where 


- ji 
2 P,f 18 


t[\R/H 2 Oj T=Tav 
The expression used for is (ref. 5) 


(56) 


(57) 


where 


N Nu ~ 


N Nu,c( q / q o) 

t , %u.c(V%) 

N Nu,FM 


(58) 


(59) 
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N Nu,C = 2 + °* 6 '| N Re 



and 


with 


q - L lnfl i C P’ f AT ) 

(60) 

q o c p,f AT ' L ' 

f Q (u - v) 2 

AT = + T - T w 

2c n w 

(61) 


N- 


Nu,FM 


iMS 1 

>-*■ 

+ 

“MS 
1 

p(u - v)2rc p>V)W N st _ 

f r (u - v) 2 _ 

2 + c p( T T w) 

r 

“!ms 

l-l. 

1 


*r.o<" - v)2 , 5 \ 

2 + c pl L w) 


(62) 


fr = 


r — 

15 

V — A 


Zn r i 


i=l 

15 

15 

v — X 

si 

n r i + Z n 

i =1 

1=! 


N 


St 


s = 


[2 + 0.7 exp(-0.7075)] 


= °- 1041 + 0.125[l - exp(-1.166S 0,406 ) 

Sl.14 L L 


(63) 


u - v 


(64) 

(65) 


15 


| 2RT ^ 


TYPICAL EXAMPLE 


Sample calculations were made for water injection into a stream tube flowing over a 
reentry vehicle at an altitude of 47.55 km (156 000 ft) with an entry velocity of 5.4 km/sec 
(17 717 ft/sec). The example given herein is intended to illustrate the use of the com- 
puter program. Necessary inputs for the program are listed in the appendix. Each 
printed output includes the inputs used. 

Figure 1(a) shows the pressure distribution used in the present example (based on 
ref. 28). With one exception, this pressure distribution was used for all the examples in 
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the figures. In order to investigate the effect of a change in this pressure distribution, 
one calculation was made by using the variation illustrated by curve B in figure 1(b) which 
was obtained from data of references 17 and 29. It should be noted that when curve B 
was used, the only input quantities changed were pressure and density, and momentum 
was not conserved between a station just upstream of injection and one just downstream 
of injection. 

The four test cases with assumed conditions are as follows: 

(1) Water and air mixture with finite rate chemistry 

(2) Air with finite rate chemistry (w* = 0) 

(3) Water and air mixture without any neutral-gas chemical reactions (frozen) 

(4) Water and air mixture with full gas chemistry but with F e = 0 (no recombina- 

tion of ions and electrons on droplet surfaces) 

The mixture temperature, density, and velocity along the body streamline are shown 
as a function of distance from the injection site in figures 2 to 4, respectively. Results 
for the first three test cases are shown for each parameter (case 4 curves are identical 
to case 1 curves) along with results obtained by using curve B of figure 1(b) (only the 
water plus air finite-rate curves are used to illustrate this effect). Figure 5 shows the 
fraction of water evaporated for both frozen and finite-rate conditions (case 4 curve is 
identical to case 1 curve). 

Figure 6 shows the reduction in electron concentration along the body streamline. 
Results for all four test cases are shown as well as the results obtained by using curve B 
of figure 1(b). The results indicate that the recombination of electrons and ions on the 
surface of water droplets is effective (for the conditions cited) in reducing electron con- 
centration in the plasma sheath surrounding a reentry vehicle. The finite-rate curve 
obtained with the curve B pressure in figure 1(b) shows a larger reduction in N e than 
that obtained with the curve A pressure of figure 1(b). 

The variation of the neutral species with and without water addition is presented in 
figure 7. The curves represent the result of finite-rate calculation which does not include 
the effects of injection-induced compression. After water injection the N and O atom con- 
centrations rapidly decrease. The importance of the removal of these atoms is that a 
condition of dynamic balance under local conditions for reaction 15 (N + O ~ NO + + e”) is 
destroyed. (See ref. 30.) It should be noted that the group of reactions (7, 8, 9, 12, 13, 

14) responsible for the increased speed of removal of the O and N atoms all contain hydro- 
gen, which was introduced by the water injection. 
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CONCLUDING REMARKS 


An analysis is presented for use in the determination of the electron concentration 
in a flowing plasma in the presence of accelerating and evaporating water droplets. For 
given input gas-flow properties and gas composition, the resulting computer program is 
used to calculate gas density, gas temperature, gas and droplet velocity, enthalpy, frac- 
tion of droplet evaporated, Mach number, and gas composition. Sample calculations indi- 
cate that a reduction in electron concentration can be achieved by injecting water into the 
plasma sheath surrounding a reentry vehicle. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., July 10, 1970. 
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APPENDIX 


DIGITAL COMPUTER PROGRAM FOR CALCULATING ELECTRON 
CONCENTRATION IN PRESENCE OF ACCELERATING AND 
EVAPORATING WATER DROPLETS 

By Barbara L. Weigel 

This appendix describes the digital computer program (D2331) developed in support 
of a stream-tube analysis used to calculate electron concentration in the presence of 
accelerating and evaporating water droplets. The program was written in FORTRAN IV 
language for the Control Data 6600 computer system with the SCOPE 3 operating system 
and library tape. 


Input 

The data input consists of one card for case identification, FORMAT (8A10), fol- 
lowed by the numeric input using FORTRAN IV NAMELIST. The NAMELIST input sym- 
bols are as follows: 


FORTRAN 

symbol 

Symbol 
in text 

Description 

$NAM1 



RHOL 

Pi 

Water density, g/cm^ 

HL 

h l 

Specific enthalpy of water, ergs/g 

DN 

d N 

Reference length given as nose diameter, cm 

EL 

L 

Latent heat of water, ergs/g 

ENTHAL 

h 2 

Initial specific enthalpy, ergs/g 

TW 

T w 

Droplet surface temperature, °K 

FRO 

f r,0 

Droplet recovery factor for zero mass transfer 
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FORTRAN 

symbol 

Symbol 
in text 

Description 

CPVW 

c P,v,w 

Specific heat of vapor at droplet surface, ergs/g-°K 

WASK 

w* 

Initial ratio of water mass flow to gas mass flow 

R 

r 2 

Initial droplet radius, cm 

HASK 

h* 

Total specific enthalpy of shocked gas, ergs/g 

U2 

u 2 

Initial gas velocity, cm/sec 

P2 

P2 

Initial gas pressure, dynes/cm2 

RH02 

P2 

Initial gas density, g/cm^ 

T 

t 2 

Initial gas temperature, °K 

VAR(l) 

X 

Distance along streamline path from injection site 
divided by djj, (s - S 2 )/d^; therefore, 

= 0 initially 

VAR(2) 

f2 

(Droplet radius divided by r 2 ) 

VAR(3) 

V 

Droplet velocity, cm/sec 

VAR(4) 

n 

Specific concentration of H^O, mol/g; y^ = 0 
at injection site 

VAR(5) 

?2 

Specific concentration of C> 2 , mol/g 

VAR(6) 

^3 

Specific concentration of N£, mol/g 


Etc. for y. where i = 4 to 14; r4,r 6 ’ y 9’ r 14 = 0 
at injection site 

x2 

DER(l) r- Derivative of r^ 
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FORTRAN 

symbol 

DER(2) 

DER(3) 

FE 

NTLUP 

VARIX 

VARDP 

PRDEL 

XTERM 

HEPS 

Cl 

CIMAX 

ELE1 


APPENDIX 


ii^text* Description 

v Derivative of droplet velocity 

y i Derivative of specific concentration of H 2 O 

Etc. for y ^ where i = 2 to 14 

F e Efficiency factor for capture of electrons by water 

droplet or catalytic efficiency 

Number of points in VARIX and VARDP arrays, ^40 

Array of x values 

Array of pressure ratios p/pg a t each x of 
VARIX array 

Print answers every PRDEL increment in x; 
PRDEL = 0.1 unless input otherwise 

Value of x at which run is to be terminated 

Percent by which enthalpy will be allowed to 
vary in two successive steps when 
-500000.0 < Enthalpy < 50000.0; HEPS = 0.1 
unless input otherwise 

Initial computing interval 0.0001220703125 unless 
input otherwise 

Maximum computing interval 0.1 unless input 
otherwise 

An array of 16 values used by the integration scheme 
to control the size of the computing interval, 

0.0 < ELE1 ^ 64.0; initially 0.5 unless input 
otherwise 
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FORTRAN 

symbol 


Symbol 
in text 


APPENDIX 


Description 


ELE2 An array of 16 values used by the integration scheme 

to control the size of the computing interval, 

ELE2 < ELE1; initially 0.1 unless input otherwise 

IBUG An array of 20 integers which must be input as 1 to 

trigger debug printouts at critical points (see 
Program Listing) 

The input subroutines are as follows: 

1. ITR1 A Newton- Raphson iteration method (see ref. 31, p. 192) 

2. CALINT A modified Runge-Kutta numerical integration scheme designed 

for chemistry problems by Charles E. Treanor (ref. 32) 

5 . 2 . . . 

3. BASIC Evaluates derivatives: r , v, y^ through y^ 

4. CHECK Controls the size of the computing interval 

5. FOFX A function subroutine to evaluate T 

6. ERROR Checks input to CALINT 

7. FTLUP A table lookup (see ref. 33, pp. 87-91) 

Output 

The output consists of the NAMELIST input, case identification, injection param- 
eters w*, r2, dj<f, and F e , initial values for p, p, h, T, u, and v, and a list of 
the species 1 to 15. A typical output is presented as a function of distance from the 
water injection site. 
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I 


(so 


tNAMl 

RHOL = 0.1E+01, 

HL = -0.1515E+12, 

ON = 0.2032E+02, 

EL = 0.245E+11, 

ENTHAL = 0.8028C264E+11, 

TW = 0.261E+03, 

FRO = 0.85E+00, 

CPVW = 0.186E+08, 

MASK = 0.35E+01, 

R = 0.78E-C3, 

HASK = 0.12624843E+12, 

IBUG = 0 1 Of Of Of Of Of Ot Of Of Of Of Of Of Of Of Of Of 

Of 0, 0, 

U2 = 0.30321E+06, 

P2 = 0.17326E-t-05 , 

RH02 = 0.17008E-05, 

T = 0 f 287E + 04 , 

PROEL = 0.1E+00, 

VAR = 0.0, 0.1D+01, 0.0, 0.0, 0.9783516D-05, 0 . 2625592 ID-01 , 0.0, 

0.58619980-04, 0.0, 0 . 20 17667D-02 , 0 . 1 43393 19D-0 1 , 0.0, 

0.101648560-05, 0 . 4929 86 5D-1 3 , 0 . 2 5 8991 78D-1 2 , 0 . 1 5938 89 ID-1 0, 

0.0, 0.1 016805 20- C5 , 0.0, 0.0, 


DER 

o o 
o o 

o o 
o o 

o o 
o o 

o o 
o o 

o o 
o o 

II 

0 • C * 0 . 0 1 
0 • C » 0 . 0i 

0.0, 0 
i 0.0, 0 

o o 

o 

. 

o 

o 

o 

EL E 1 

= 0 .5E+0C * 0 • 5E + 00* 0.5E+OC, 

0.5E+0G, 0.5E+00, 0.5E+0C, 

0 . 5E+0G * 0.5E+00, 0.5E+OC, 

0 . 5E + CC , 
0.5E+00* 
0.5E+00, 

0.5E+00, 

0.5E+00, 

0.5E+00, 

0.5E+00, 

0.5E+00, 

0.5E+00, 

0.5E+00 

0.5E+00 

ELE2 

= 0 . 1 E +00 , 0.1E+00, 0.1E+0C, 

0 • IE + 00 » 

0.1E+00, 

0.1E+00, 

0.1E+00 
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O.IE+OO, O.IE+OO, O.IE+OO, O.IE+OO, O.IE+OO, O.IE+OO, O.IE+OO, 

O.IE+OC, O.IE + OO, O.IE+OC, O.IE+OO, O.IE+OO, O.IE+OO, 

Cl = 0.1220703125E-03, 

CIMAX = O.IE+OC, 

XTERM = 0.45E+C1, 

HEPS = O.IE+OO, 

FE = 0.3E-01, 

NT LUP = 27, 

VARIX = O.O, 0.1018E + 00, 0.2011E+00, 0.3026E+00, 0.4003E+00, 

0.5058E+00, 0.611 2E+00 , 0.703E+00, 0.8046E+00, 0.9062E+00, 

0.1 00586+01 , 0.11034E+01, 0.13066E+01, 0.15097E+01, 0.17128E+01, 

0.19003E+01, 0 .2 1034E+01 , 0.23065E+01, 0.25097E+01, 0. 271286+01, 

0.29003E+01, 0.31 034E+01 , 0.33065E+01, 0.35097E+01, 0.37128E+01, 

0.39003E+01, 0.4D878E+01, I, I, I, I, I, I, I, I, I, I, I, I, I, 

VARDP = 0.1E+01 , 0.923E+00, 0.858E+C0, C.8E+00, 0.75E+00, 0.705E+00, 

0.667E+00, 0.639E+00, 0.612E+00, 0.591E+00, 0.573E+00, 

0.558E+00, 0.536E+00, 0.522E+00, 0.512E+00, 0.5C5E+C0, 

0.498E+00, 0.489E+00, 0.478E+C0, 0.465E+00, 0.451E+C0, 

0 . 436E+00, 0.42E+00, 0.4C8E+00, 0.401E+00, 0.401E+00, 

0.401E+00 , I, I, I, I, I, I, 1, I, I, 1, I, I, I, 

SEND 

D2331 - STREAMTUBE//REPORT // PROGRAM 
TEST CASE FOR REPORT 


R= 7.80000000E-04 DN= 2 .03200000E+01 

RHO= 1 • 70080000E-06 ENTHALPY^ 8.02802640E+10 

U2= 3.0321C000E+05 V2= 3.45719769E+02 



to 

U1 


APPENDIX 


CO 

Oi 


0 + ' 

13 

H + 

14 

E- 

15 


X= DISTANCE ALONG STREAMLINE FROM INJECTION SITE/NOSE DIAMETER 
RHO=G AS DENSITY, G/CM **3 

R HAT-DROPLET RADIUS DIVIDED BY INITIAL DROPLET RADIUS 

LAMBDA SUB D=DEB YE LENGTH , CHARACTER I ST IC DISTANCE OF FIELD BEYCND WHICH THE EFFECT OF A CHARGE IS NOT FELT, CM 
TE MPERATURE=GAS TEMPERATURE, DEG. K U=GAS VELOCI TY , CM/SEC 

Z- COM PRESS IB IL I TY FACTOR M BAR=ME AN MOLECULAR WEIGHT FOR IONS,G/MOLE 

S=DI ST ANCE ALONG STREAMLINE FROM INJECTION SITE, CM T IME=EL APSED TIME AFTER WATER I NJ ECT ION , S EC 


V=DROPLET VELOCITY, CM/SEC 
A SUM=SUM OF MASS FRACTIONS OF SPECIES 
PRESS URE=GAS PRESSURE , DYNES/CM **2 


ENTHAL PY=SPECI FIC ENTHALPY OF GAS , ERG/ G 

W SUM=$UM OF ALL THE GROSS REACTION RATES, MOLE/CM** 3 -SEC»=C.O 
MACH NUMBER 


N SUB E=ELECTRON C 0 NC 6 NT RAT I ON , EL ECT RONS/CM * *3 

FE= EFFICIENCY FACTOR FOR CAPTURE OF ELECTRONS BY WATER DROPLET OR CATALYTIC EFFICIENCY 


F=FRACT I CN OF WATER EVAPORATED 

CAPITAL U=NQRMA L I ZED RELATIVE VELOCITY BETWEEN DROPLET AND GAS 


RHO 

TEMPERATURE 

U 

V 

ENTHALPY 

N SUB E 

F 

R HAT 

Z 

M BAR 

A SUM 

W SUM 

FE 

CAPITAL U 

LAMBDA SUB D 

s 

TIME 

PRESSURE 

EM 



1 • 69794 50 3 E -0 6 

2 . 87 442 936 E +03 

3 . 032 10000 E +05 

3 . 45719769 E +02 

8 . 02800688 E +10 

1 . 03986930 E +12 

0 . 

1 .OOOCOOOOE+OO 

1 . 23098765 E+ 0 C 

2 • 999 9 77 6 3 E + 01 

9 . 94 94 439 7 E - 01 

- 3 . 87771583 E -18 

3 • OOOOOOOOE -02 

2 • 12379678 E +00 

3 . 62773481 E - 04 

0 . 

0 . 

1 . 73 2 60 000 E + 04 

2 . 50752 A 29 E +00 




GAMMA, SPECIFIC CONCENTRATION OF SPEC IES , MOLES/G,GN RETURN FRGM INTEGRATION 


I 0. 

6 0 . 

U 4 . 92986 50 OE - 1 4 


9 . 78351 60 CE -06 
2 . 017667 COE-C 3 
2 . 58 99 1780 E - 13 


2 . 6255921 0 E -02 
1 . A 3393 1 9 OE - 02 
1 . 5 S 388910 E -11 


0 . 

0 . 

0 . 


5 . 861 99800 E-C 5 
1 . 01648560 E -06 
1 . 016501 85 E -06 


X, MOLE FRACTION OF 

0 . 

0 . 

1 . 15498567 E -1 2 


SPECIE 

2 .29 2 11 566 E — 04 
4 . 727059 40 E -02 
6 . 06 774818 E - 12 


6 . 15132716 E — 01 0 . 
3.3594648 0 E -01 0 . 
3 . 73421801 E -10 0 . 


1 . 37336898 E - 03 
2 . 39145730 E-C 5 
2 . 381 4953 7 E -05 


W, NET RATE OF PRODUCTION OR DISAPPEARANCE OF SPEC I E S , M 0 L£/CM** 3 -SEC 
1 0 . - 4 . 10506688 E-C 6 2 . 795 228 1 2 E -06 0 . 

6 0 . - 6 . 89624362 E-C 6 6 . 88302419 E-G 6 0 . 

11 7 . 13939703 E -1 2 - 2 . 19 758 12 IE - 1 2 2 . 13342 736 E -08 0 . 


9 . 99898901 E - 06 
- 8 . 69321371 E -06 

0 . 


LOWER CASE K SUB J 
8 . 12807257 E +06 
2 . 35421 897 E + 1 3 
4 . 30167004 E+ 0 6 
1 . 96440443 E +0 3 


, J= 1 , 45 ,REACTICN 
3 . 51 01 2965 E +06 
9 . 53710703 E +09 
7 . 31 101221 E +01 
4 . 619277 1 3 E + 1 2 


RATE CONSTANTS FOR A PARTICULAR REACT ION, CM** 3 / MOLE-SEC 

6 . 84643620 E -02 1 . 06394682 E +08 1 . 965 55567 E +04 6 . 47714336 E +07 
5 . 76345143 E +10 7 . 3 1 42 1473 E +07 8 . 475 6622 5 E +1 1 3 . 79601923 E + 12 
8 . 45881519 E -13 1 . C 6378657 E -01 5 . 63 76 1 1 05 E -1 1 2 . 40840919 E -09 
1.0 3502394 E +1 1 3 . 22759684 E +1 1 1 . 61271 5 57 E +05 6 . 58923026 E +06 


2 • 40555 199 E +1 2 
8 . 72018 00 OE+IO 
2 . 50886332 E - 09 
2 . 251 52743 E +07 
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8 . 052 19631E+1 2 

2.18779306EM3 

9.30829116E+14 

6 • 10644081E+02 

3.622979776+10 

3.73281198E+12 

1.31103829E+D9 

1.01CA6708E<-02 

2.20AA136CE+10 

4 • 46320053E+09 

1 • 27294680E+07 

9.24264027E+02 

2. 95285058E+1 2 

2.95285058E+1 2 

2 .95285058E+1 2 

2.95285C58E+12 

2. 95285C58E+1 2 





UPPER CASE K SUB J, 

J=l. AO, EQUILIBRIUM CONST ANTS , MOL ES/C M**3 OR DIMENSIONLESS 



1.52911 A82E-0 8 

2 • 20415997E-08 

1. 486664536-16 

A.63367A23E-C8 

1 .70046476E-11 

2.83 528 1 81E-08 

7.77A0A193E-01 

1. 63A2905AE+00 

5,99 751 5846— 04 

7 * 714 79739E-04 

8.7A2 69537E-06 

5 . 3931 6698E-01 

3. 30000503E-01 

1.817A2 627E-D3 

6.877 6A 1 A7E- 1 0 

9.07544680E-15 

1.34921449E-30 

1 • 1695 1869E-20 

9 .01794061 E-29 

7.65395155E-28 

9.07321365E-28 

1. 15364936E-10 

1. 49614479 E-02 

1 . 76276853E-03 

1.A87030AAE-03 

7. 7108135 7E- 09 

6. 5AA53118E-08 

7.75807A93E-08 

1.178207 17E-0 1 

9. 93 90 8 1086- 02 

8.43576691E-01 

5.9A8736AAE-12 

5 . 84348 2 1 3E-04 

7.A8571339E-03 

1 • 31955800E-05 

A. 62A57265E-1 2 

9 *089629626-05 

5 • 96C97962E-0 5 

1.C300706AE-06 

1.40138352E-11 



CONTRIBUTION OF INDIVIDUAL REACTIONS 

TO GROSS REACTION, W ,MOLE-CM**3/G**2 -SEC 



AA 

AB 

AC 

AD 

AE 

AF 

AG 

AH 

A I 

AJ 

AK 

AL 

AM 

AN 

AO 

AP 

AQ 

AR 

AS 

AT 

AU 

AV 

AW 

AX 

AY 

AZ 

BA 

BB 

BC 

BO 

BE 

BF 

BG 

BH 

BI 

BJ 

8K 

BL 

BM 

BN 

BO 

BP 

80 

BR 

BS 

BT 




0. 

2.37165123E+03 

1. 3587275AE+02 

0. 

2. 423 6448 1E+C3 

0. 

0. 

0. 

0. 

1 . A262503 3E + 06 

9.619664796+05 

0. 

0. 

0. 

6.33813958E+03 

1.0606A516E-0A 

5 . A225A738E-1 A 

1.62214013E-05 

2.84062243E-13 

8 . 7991 8532E- 1 1 

0. 

-3.21929658E-06 

-7. 01037076 E—CA - 

-1.8080Q066E-03 

0. 

-1.322332 59 E- 05 

-1.97081906E-03 

0. 

-5.1AAA9010E-03 

0. 

0. 

2 • 51201831E-04 

9. 15684008E-04 

-7.AA7050A5E+03 

-2.61861255E+00 

0. 

5.9391175AE-0A 

0. 

1 .82291010 E-04 

0. 

3 . 00153009E+06 

1 . A55715A7E-0 1 

7.6A76A028E-01 

A. 70651 636E+C1 

0. 

0. 




GAMMA CORRECTED TO 

MAKE ASUM=1.0 






0. 

9.7 8 351 60 CE- 06 

2.62559210E-02 

0. 

5.86199800E-C5 



0. 

2.01 7667C0E-03 

1 . A33931 90E-02 

0. 

1 • 0 164 8560E-06 



A. 92986 500E-1 A 

2.5899 1780E- 13 

1. 5938891 OE-11 

0. 

1.016 501 8 5E-06 



ASUM= 9.9A9AA397E 

-01 






FGAM 1 1 1 = 0. 







GAMMA, COMPOSITION AFTER WATER ADDITION AND CORRECTION 





0. 

9. 7 8 35 1 60 OE- 06 

2 . 625592 1 9E-02 

0. 

5 • 861 998 00 E-C 5 



0. 

2 . 0 1 7667 COE- 03 

1.433931 906-02 

0. 

1.0 1 648560E-06 



4 . 92986 500E- 1 4 

2.5899 17806-13 

1. 593889106-11 

0. 

1.0 16 501 85E-C6 



ASUM= 9.94944397E 

-01 






RHO 

TEMPERATURE 

U 

V 

ENTHALPY 

N SUB E 

F 

R HAT 

Z 

M BAR 

A SUM 

W SUM 

FE 

CAPITAL U 

LAMBDA SUB D 

S 

TIME 

PRESSURE 

EM 



1.72137253E-06 

2 • 40 98 532 2 E+03 

1. 61103489E+05 

1 . 819 1 8868 E + 04 

9.57785969E+09 

2 • 78A2 1A66E+1 1 

1.69A2A865E-01 

9 . 391 28313E-0 1 

1.09 199 131 E+00 

2 .999999956+01 

9. 9682 6338 E- 01 

8 . 972 4008 5E-1 7 

3 . OOOOOOOOE-02 

1 • 09 A A9 68 9 E+00 

6.419379166-04 

2.06 87 1094E+00 

1 • 85793400E-04 

1 .5991322764-04 

1.30061 27 7E +00 




GAMMA, SPECIFIC CONCENTRATION OF SPECIES ,MOLES/GtON RETURN FROM INTEGRATION 


tS3 


APPENDIX 


CO 

00 


1 1.95926887E-02 

4 . 87 344039E-04 

1.65312253E-02 

1 .69369741E-04 

1 . 2031 1686E-03 



6 2.28948921E-03 

1.53 114292E-07 

5. 61645665E-03 

4. 61557202 E-04 

2.68460796E-07 



11 9 • 83179955E- 1 B 

2.22 850465E-20 

8. 048092 49E-15 

7. 44084463 E-16 

2.68460805E-07 



X.MOLE FRACTION OF 

SPECIE 






4.22694058E-01 

1.05 1 39 9 49 E— 02 

3. 5664583 IE— 0 1 

3.65399483E-03 

2 . 59 5612 87E-02 



4.93936030E-02 

3. 30329861 E-06 

1.21169835E-01 

9.95766789E-03 

5.791 79230E-06 



2-12111 942E-1 6 

4, 807 79 177 E— 19 

1.73630119E— 13 

1 .6052931 IE-14 

5.79179249E-06 



W » NET RATE OF PRODUCTION OR DISAPPEARANCE OF SPECIES* 

M0LE/CM**3-SEC 




1-1-66511595 E- 04 

1. 00 79 2367 E- 04 

2. 69013807E-09 

2.95296302E-05 

3 • 1 5594067E-08 



6 2 - 02680963 E- 04 

-2. 330 2098 8 E- OB 

-2 . 3777202 5E-04 

7.12829662E-05 

-1. 3637584 0E-08 



11-4. 01122446 E-l 8 

-8.13721508E-21 

-1.39217644E-15 

1 .04353378E-16 

0 . 



LOWER CASE K SUB J, 

J=l, 45, REACTION 

RATE CONSTANTS FOR 

A PARTICULAR REACT ION , CM**3/ MOLE- 

SEC 


1. 95426066 E+05 

7.7651 3317E+04 

4.09996994E-05 

3. 85676599E+06 

1 .49340249E+02 

2.53 768877EFQ6 

1.2381 8433E+1 2 

1.99415828EH3 

1.68996435E+09 

1 .81673860E+10 

5.25419083E+06 

4 . 39668816E+1 1 

2. 12856604E + 12 

2.52317334E+10 

4 . 56004 330E+0 5 

7.21 639494E— 0 1 

4. 19 1 561 60E-1 8 

7 . 27298 171 E-05 

6. 277 03 175 E-16 

5. 53027187E-14 

5.80650034E-14 

1.303 64 309EF01 

1. 38248889E+ 12 

2.04520C52E+10 

6.32836138E+10 

2.40479731E+03 

2 • 02634026E+05 

6.9774481 4EF05 

3.57487553E+12 

9. 63765 43 2 E+ 12 

8 . 456 8371 6E+ 14 

3 . 13449339E+0 0 

5. 8474311 7E +09 

1 . 4631 7120E+12 

1 . 10898205F + 08 

2.91576460E-01 

3.24 057932E + 09 

3. 83458590E+08 

4.43395622E+05 

4.35530839E+00 

1. 51 149415E+10 

1.51149415E+10 

1.511 49 415E+10 

1.51 149415E + 10 

1*51 14941 5E+1 0 





UPPER CASE K SUB J, 

- J=l»40 f EQUILIBRIUM CONST ANT S r MOLES/C M**3 OR DIMENSIONLESS 



2 . 792 76 7 84E- 1 0 

4. 2042935 8E-1C 

7.53940473E-2O 

1.39064380E-09 

1 .08719305E-13 

9.02462977E-10 

4.6 5868816E-01 

1.54094277EV00 

1 . 20469546E-04 

2.5859113 5E-04 

6.93481689E-07 

3.09460655E-01 

2.00 82 5 534E-01 

5.99871 658E-04 

6 . 03392 29 7E- 1 1 

7.22874233E— 17 

5. 4500992 5E -36 

6. 560039 12E-24 

8.396 11722 E- 34 

1 .48248981E-32 

1.75885134E-32 

8. 30802 857E-1 3 

6.49121386E-03 

3. 67631481 E— 04 

3.09866964E-04 

1 . 279 888 28E-10 

2.25987952E-09 

2.6R115983E-09 

5.66352441E-02 

4, 77363665 E- 02 

8.42873855E-01 

3.30967765E-14 

1.068338 99 E-04 

3. 25 87 4433E-03 

1.19801701E-06 

1. 54187561 E- 14 

1. 46 45 3721 E- 05 

5.75077 841 E-06 

3 .92755047E-08 

7. 6776870 7E -14 



CONTRIBUTION OF INDIVIDUAL REACTIONS TO GROSS REAC TION , W ,M0LE-CM**3/G 

**2 -SEC 



AA 

AB 

AC 

AD 

AE 

AF 

AG 

AH 

A I 

AJ 

AK 

AL 

AM 

AN 

AO 

AP 

AQ 

AR 

AS 

AT 

AU 

AV 

AW 

AX 

AY 

AZ 

BA 

BB 

BC 

BD 

BE 

BF 

BG 

BH 

81 

BJ 

BK 

BL 

BM 

BN 

BO 

BP 

BQ 

BR 

BS 

BT 




-1. 18477613E+02 

4.63 107277E+02 

9. 85799910E-07 

1 . 68632752E+01 

8.5923887 8E-C2 

3. 123151 80E+02 

3.38976000E+07 

-5.29427412E+06 

3. 97916 64 OE + 03 

-1. 17519236E+05 

9.07872035E+02 

-4.0934534 4E +07 

-1.51389271E+07 

1 • 21040810E+05 

5.44666487E+02 

2.63494269E-08 

9. 380 82 62 9E— 19 

3 • 69252596E-07 

2.066931 17E-21 

3.72464048E-15 

3.04750336E-16 

1. 27754582 E- 07 

—2 .56935892E— 10 

3.50948151E-07 

1 .48347306E-07 

4.049! 5! 08E-10 

5. 62684883E-04 

1 .465I4557E-04 

3.49517822E-09 

9.78621642E-10 

1 . 92826 32 3E-04 

-3.69762127E-09 

1 • 144603 1 6E-09 

-4.06 997830E-04 

5.53109961E-07 

9.28710681E-10 

— 2.89 747700E-08 

3.05879576E-09 

2 • 43 69 3 563 E- 10 

2 .4570261 4E -08 

4.05776924E+03 

1 • 48607075E-07 

3 • 36837 174E- 1 0 

1.21 646 44 7 E- 04 

1. 12467932E-C5 

0 . 




GAMMA CORRECTED TO 

MAKE ASUM=1 . 0 






1 .9577B953E-02 

4.88 247 4 5 5 E~ 04 

1. 65323972E-02 

1 .65707045E-04 

1.2043461 5E-03 



2.29057128E-03 

1.53 1 25 1 45E-07 

5. 626 8681 8 E— 03 

4. 51575824 E-04 

2. 6873509 6E -07 




APPENDIX 



9.83249648E-18 2. 22866262E-20 8 . 06301 1 67E-1 5 7. 2T993308E-16 2 . 684608 05E-07 

ASUM = 9.96826337E-01 

FGAMI 1 }= 2.002 505 18E-0 2 

GAMMA, COMPOSITION AFTER WATER ADDITION AND CORRECTION 

1.97973453E-02 A. 84000003E-04 1 . 64487353E-02 1 . 6028 1028E-04 1 . 19591424E-03 

2 . 26740488E-0 3 1 . 52 3502 59E-07 5 . 57791789E-03 4.36789138E-04 2 . 66 85362 0E-07 

9.78273935E-1 8 2.21738452E-20 7.99286844E-15 7.04155432E- 16 2 . 66 853629E-07 

ASW= 9. 96 8 42 2 92 E -01 

X RHO TEMPERATURE U V ENTHALPY N SUB E F 

R HAT Z M BAR A SUM W SUM FE CAPITAL U 

LAMBDA SUB D S TIME PRESSURE EM 

.201073 1 .72547798E-06 2.22606350E+03 1 . 43560973E+05 2 .32739587E+04 -3. 70713870E+09 1 . 87574931EH 1 2. 14887789E-01 

9 • 22256962E-0 1 1 . 06640393E+00 2 . 99999999E+01 9 .971 14558 E-01 4. 01 790934E-17 3 . 00000000E-02 9. 58500428E-01 

7 . 51675626E- 04 4.08581044E+00 2. 82315746E-04 1 .48660103E+04 1 . 14438 1 13E+00 

GAMMA, SPEC IFIC CONCENTRATION OF SPEC I ES , MOLE S/G ,0N RETURN FROM INTEGRATION 
1 2 . 2282896 5E-0 2 1 . 10 026245E-03 1. 50313788E-C2 1 .64965814E-04 1 .09305398E-03 

6 2 . 66234168E-03 4.67983605E-C8 3 . 82319735E-03 3 . 7683 G759E-04 1 .80434014E-07 

11 8. 16984892E-19 9.42418774E-22 1 . 18530402E-15 1 . 295621 71 E- 1 6 1 . 80434016E-07 

X, MOLE FRACTION OF SPECIE 

4.788382 1 OE-O 1 2. 3643591 5E-02 3 . 2301 0005E -01 3 . 54495813E-03 2 . 34886883E-02 

5.72111854E-02 1.00565217E-06 8. 21568673E-02 8 . 09773388 E-03 3 . 87735501 E-06 

1 .7556226TE-17 2 .02516813E-20 2 . 54710537E-14 2 .78416757E-15 3.87735504E-06 

W.NET RATE OF PRODUCTION OR DISAPPEARANCE OF SP ECI ES ,M0LE/CM**3- SEC 
1—7. 80187245 E- 05 7.420421 78E-05 8. 67359525E-10 2.87519424E-05 7.88704862E-09 

6 5 . 307907 53E-0 5 -4. 06 1 04548 E-09 - 1 . 23471 1 1 3E -C4 4.54544889E-05 -5 . 560722 1 9E-09 

11 — 1.43 26962 9E-19 -2. 14307227E-22 -1 .49570299E-16 1 . 23 21 91 33 E-l 8 0. 

LOWER CASE K SUB J , J=1 ,4 5, REACTI CN RATE CONSTANTS FOR A PARTICULAR REACT ION, CM**3/M0LE-SEC 

2.86875975E+04 1 . 09 713 539E+04 9. 16271548E-07 7. 0135 1401 E+05 1.22217840E+01 4. 801 74630E+05 8.86695751E+1 l 

1 • 84188940E+ 13 7.01932702EF08 1 .0L2732C0E+10 1 . 37598847E+06 3. 161 12995E+11 1 . 59247775E+12 1 . 34633227EF1 0 

1.45681 828E*05 6 . 8 56 36 846E-02 8. 218 12254E-2 1 1 .76678632E-06 1 . 861 041 70E-18 2. 37318709E-16 2 . 50177146E-16 

1 .01124739E+00 1 . 1965122 2E+ 12 8.98108419E+09 2 . 76796661 E+ 10 2. 820 86437E+02 3.44032773E+04 1 . 18929913E+05 

2.37379145E+12 6.37419902E+12 8. C9569402E+14 2 .1324T248E-01 2.31559622E+09 9. 11689839E+11 3. 15709573E+07 

1.47803385E-02 1 . 22347779E+09 1 . 10036415E+08 8.02262767E+04 2. 83631413E-01 8 . 87205956E+09 8.87205956E+09 

8.87205956E+09 8. 87205956E+09 8 . 87205956E+09 

UPPER CASE K SUB J , J=1 , 40, EOUI L I BR IUM CONSTANTS, M0LES/CM**3 OR DIMENSIONLESS 

3 . 59375628E- 1 1 5. 52 524322E-1 1 1 . 55631902E-21 2 . 31 104058E- 10 8 . 19001 5 52E- 1 5 1 . 54591091E- 10 3. 57410197E-01 

L.49493776E+00 5. 29785739E-05 1 .48229050E-04 1 .900 26383E-07 2.3246B525E-01 1. 55503815E-01 3.40689866E-04 

1.77701694E-11 6. 28 101053E- 18 9. 77525617E-39 1 .45537963E-25 2 . 29731 047E-36 5. 89288371E-35 6.99946088E-35 

6. 71663664 E- 14 4 . 25 5 087 1 2E-03 1 . 65882387E-04 1 . 39657273E-04 1 . 57849568E-11 4. 04903544E-10 4. 809371 1 9E-1 0 

3.89844867E-02 3 . 28 2 1 248 8E-02 8 . 41905 372E-01 2 .33978915E- 15 4.46586350E-05 2 . 13077541E-03 3. 53458112E-07 

8.36264502E-1 6 5.77863342E-06 1 . 73882295E-06 7.40808099E-09 5.37777482E-15 

CONTRIBUTION OF INDIVIDUAL REACTIONS TO GROSS REACTION, W,M0LE-CM**3/G**2-SEC 


to 

co 




i 
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CO 

o 


AA 

A0 

AC 

AD 

AE 

AF 

AG 

AH 

A I 

AJ 

AK 

AL 

AM 

AN 

AO 

AP 

AQ 

AR 

AS 

AT 

AU 

AV 

AM 

AX 

AY 

AZ 

BA 

BB 

BC 

BD 

BE 

BF 

BG 

BH 

81 

BJ 

BK 

BL 

BM 

BN 

30 

BP 

BQ 

BR 

BS 

BT 




3.45581001E+01 

2.32491161E+02 

1.02892153E-07 

2.92188523E+01 

2. 08170831E-02 

2.99828655EF02 

2.48845051E+07 

7 . 44175024E+05 

1.36165816E+03 

-3.83037382E+04 

2.91326639E+02 

-1. 72918 51 8E +07 

-8.874150! 3EF06 

3.87816483EF04 

2 .66901 244E+0 2 

1 • 6091 5184E-C9 

3.85837771E-20 

1 .23047387E-07 

4. 2 88 709 1 3 E— 23 

2.68 151427E-16 

2 . 60 140224E- 1 7 

1.07C23351E-08 

-6.19850588E-12 

9.09621298E-09 

7.11 1 8 6044 E- 09 

1.60268095E-11 

8.63502549E-05 

2 .69341 111 E-05 

8* 77174149 E-l 1 

3.03 2 14270E- 1 1 

2 • 84901 746E-05 - 

-5 • 26 02 09 12 E- 11 

4. 61797504E-11 

-7.51113626E-05 

1. 2405551 4E -08 

4.33404250E-11 

-1.36560064E-09 

1.43 194251 E-l 0 

9. 72631 86 IE- 12 

1 . 099 1 l 765E-09 

1.60082132E+03 

7. 24833862 E-09 

8.36119549E-1 2 

1.05160878E-05 

1. 14948330E-06 

0. 




GAMMA CGRRECTED TO 

MAKE ASUM=1 *0 






2. 22796429 E- 02 

1 • 10 076884E-03 

1 . 50314954E-02 

1 • 64141629E-04 

1 .09 3 3 26 2 4 E-03 



2.662712 50E-03 

4. 679872 37 E- 08 

3 • 8249569 5E— 0 3 

3.7494 8076E-04 

1 .80478958E-07 



8.16991231E-19 

9 .42 42 608 6E— 22 

1. 18584955E-15 

1 .28914865E-16 

1.80434016E-C7 



ASUM= 9, 971145 58E 

-01 






FGAM ( 1 )= 2.241 17214E-02 






GAMMA, COMPOSITION AFTER HATER ADDITION AND CORRECTION 





2.23376245E-02 

1 • 09 783970E-03 

1 . 501 0999 9 E -02 

1 .62751753E-04 

1.09107891E-03 



2. 65471755 E- 03 

4. 67 3491 3 3E- 08 

3. 81477875E-03 

3 . 717731 88E-04 

1. 80107983 E-07 



8.15877258E-19 

9.411410 83E- 22 

1. 1 8269400E-1 5 

1 .27823274E-16 

1 .80107984E-07 



ASUM= 9.97118468E 

-01 






X RHO 

TEMPERATURE 

U 

V 

ENTHALPY 

N SUB E 

F 

R HAT 

Z 

M BAR 

A SUM 

W SUM 

FE 

CAPITAL U 

LAMBDA SUB D 

S 

T IME 

PRESSURE 

EM 



.302636 1.66198967E-06 

2 . 14402 1 18E+03 

1 . 35374 84 2E +05 

2 . 67903644E+04 

-1.14768154E+10 

1. 33 692994E + 11 

2. 4439810 8E-01 

9.10419934E-01 

1 . 05570381 E + OQ 

2 • 99999999E+0 1 

9.97275186E-01 

7. 9 49 77468 E- 18 

3. OOOOOOOOE-D2 

8.81648663E-01 

8. 73794489 E- 04 

6.14956044E+00 

3 . 6436682 3E -04 

1.3860 482 1 E+04 

1 • 05 4 24 749 E +00 



GAMMA, SPECIFIC CONCENTRATION OF SPEC I ES ,MOLES/G,ON RETURN FROM INTEGRATION 



1 2.40430942E-02 

1.48 67413 2E-C3 

1.42C49716E-02 

2.65708441E-04 

1.0107781 IE— 03 



6 2 . 62483 52 7E-0 3 

2 • 2606679 4E-C8 

2 . 70472 592 E-03 

4.27666964E-04 

1 .33516017E-07 



11 2. 18522759E-19 

1.48 628 809 E— 22 

3 • 60464867E- 1 6 

6.17403821E-17 

1 • 335 16017E-07 



X , MOLE FRACTION CF 

SPECIE 






5. 14C83926E-01 

3 . 17 891 6 1 8E— 02 

3. 03727444E-01 

5.681 31694E-03 

2 • 161 22258E-02 



5.61236258E-02 

4. 83370835 E-C7 

5. 7831829 3E -02 

9. 14427694E-03 

2.85480885E-06 



4.67240351E-1 8 

3 • 177947 12 E- 2 1 

7 • 707377 12E-1 5 

1.32011869E-15 

2.8 548088 6E-06 



H.NET RATF OF PRODUCTION OR DISAPPEARANCE OF SPECIES, 

M0LE/CM**3-SEC 




1—4. 25930438E-0 5 

4.58174821E-05 

3 . 7263366 8E-1 0 

2 . 19186816E-05 

3. 045730 77 E- 09 



6 1.74502 734E-0 5 

-9.51025472E-10 

-6. 64923996 E-05 

2. 3898451 0E-05 

-2 • 839972 63E-Q 9 



11-2.45082712E-20 

-9.08483583E-24 

-2.30311294E-17 

6 .28855450E-19 

0 . 
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LOWER CASE K SUB J. 

J= 1,45, RE ACTION 

RATE CONSTANTS FOR 

A PARTICULAR REACT ION , CM**3/M0LE- 

-SEC 


1 • 092 13026E+0 A 

4.10 122674E+03 

1. 3 57 7672 4E -07 

2.97569119E+05 

3.47366234E+00 

2.07829410E+05 

7 • 507 68476E+1 1 

1.77192155E+13 

4.52199069E+08 

7 . 56211154E+09 

7.03308804E+05 

2.601 88677E+ 1 1 

1.37 827627E + 12 

9, 834R8526E+09 

8. 2 27 17132 E+ 04 

2. 10713274E-02 

3. 6037121 2E-2 2 

2 .73905215E-07 

1.00426544E-19 

1. 54168477E-17 

1.62850135E-17 

2 . 80714 198E-0 1 

9.54236712E+11 

5.94966453E+09 

1 .83003717E+10 

9 • 636 16226E+01 

1.41 479 302E+04 

4. 90050326 E+ 04 

1 .93 5 14294E+ 1 2 

5. 18596309E+12 

7 . 9292691 6E+ 14 

5.54242679E-02 

1 .45646774E+09 

7. 19871050E + 11 

1 . 68270977 E+07 

3.31425981E-03 

7.51339005E+C8 

5 . 88831 81 7E+07 

3.40625144E+04 

7.211 86369E-02 

6.56398738E + 09 

6.56398738E+09 

6. 56398738 E+ 09 

6.56 398738E + C9 

6.56398738E+09 





UPPER CASE K SUB J. 

J=1»40,EQUILI BRIUM CONSTANT S , MOL ES/C M**3 OR DIMENSIONLESS 



1.28320452E-1 1 

1.99333935E-11 

2 . 2 1878201E-22 

9.38571703E-11 

2 .235587626-15 

6. 37492052E-11 

3. 12684581E-01 

1.47228769E+00 

3. 50684783 E— 05 

1. 1215288 7E-04 

9.92482690E-08 

2.01289493E-01 

1 . 3671 8859E-01 

2.56500666E-04 

9.66438565E-12 

1.85321588E-18 

4. 1 1188205E-40 

2. 16055809 E-26 

1.19355725E-37 

3.69198614E-36 

4.38857690E-36 

1.90315737E-14 

3.44506479 E- 03 

1. 11373171E-04 

9.3695 1123E-05 

5. 52430067E-12 

1.70881133E-10 

2.03122375E-10 

3 • 23283 242E- 0 2 

2.71969087E-02 

8.41271836E-01 

6 . 19 566266 E- 16 

2 ■ 88088251 E-05 

1.72 175429E-03 

1. 91757236 E- 07 

1 .93728 818E-1 6 

3.62571488E-06 

9 . 53754202E-07 

3. 20853022 E- 09 

1.41698680E-15 



CONTRIBUTION OF INDIVIDUAL REACTIONS TO GROSS REACTION, W,M0LE-CM*«3/G**2-SEC 



AA 

AB 

AC 

AD 

AE 

AF 

AG 

AH 

A I 

AJ 

AK 

AL 

AM 

AN 

AO 

AP 

AQ 

AR 

AS 

AT 

AU 

AV 

AW 

AX 

AY 

AZ 

BA 

BB 

BC 

BD 

BE 

BF 

BG 

BH 

81 

BJ 

BK 

BL 

BM 

BN 

BO 

BP 

BO 

BR 

BS 

BT 




6.1 9820965E+0 1 

1 • 16708864E+02 

2.421 88975E-08 

4. 1375 1478E+01 

7. 22062154E-03 

2. 67607042E+02 

1 .656 8732 8E+07 

7.75889034E+05 

5.69683802E+02 

-1.840765 43 E+ 04 

1 • 34904162 E+02 

-8.26071103E+06 

-7. 14106883E+06 

1 • 81991 129E+04 

1.51755 001ET0 2 

3.3173 85 28E- 10 

5.67415476E-21 

5.01490581E-08 

3 • 705 12889E-24 

4.45957552E-17 

6. 78779601 E-l 8 

2.72553317E-09 

-6.46315259E-13 

1. 10954004E-09 

2 • 20366770E-09 

2.32965535E-12 

2.50568285E-05 

1.22577399E-05 

8. 29403891E- 1 2 

4.88216327E-12 

1. 2888438 OE -05 - 

-2. 4759 7 87 5E- 13 

6. 73427469E-1 2 

-3. 19725253E-05 

1.07522540E-09 

6. 484926 15E-1 2 

-2.27959306E-10 

2 . 041 1 0787E-1 1 

1 .42047163E-12 

1 . 793 1 7356E-10 

8. 76397451E+02 

1.43 438 063E-09 

9.7 5597629E- 1 3 

2. 3660 86 84E— 06 

4. 05263 08 9E-07 

0. 




GAMMA CORRECTED TO 

MAKE ASUM= 1.0 






2 • 40373549E-0 2 

1. 48794069 E- 03 

1 . 42C545 15E-02 

2 • 6342 2797 E- 04 

1.011 22894E-03 



2.62550002E-03 

2.2607 4433E-0 8 

2.7069C785E-03 

4 • 23 988 142 E- 04 

1 .33575567E-07 



2 • 18530142E- 19 

1.4863383 IE- 22 

3. 60755659E-16 

6 • 12 09 2868 E- 17 

1.33516017E— 07 



ASUM= 9. 972751 85E 

-Cl 






FGAM< 1 )= 2.42293046E-02 






GAMMAtCOMPOSITION AFTER WATER ADDITION AND CORRECTION 





2. 41137262 E- 02 

1.48245061E-03 

1.41788466E-02 

2.59 88 9254 E- 04 

1 . D0 835517E-03 



2.61431084E-03 

2.25651026E— C8 

2 . 696 920 1 3E-03 

4 . 18300782 E-04 

1. 3319 59 65 E- 07 



2.1 81 20865E-1 9 

1.48 355460E-22 

3 • 59424572E-1 6 

6. C3 882279 E- 17 

1.33195965E-C7 




ASUM= 9.97280250E-01 
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The definitions of the output symbols are given as follows: 


FORTRAN Symbol 

symbol in text 


Description 


X x 

RHO p 

RIIAT r 

LAMBDASUBD Ad 
TEMPERATURE T 
Z Z 

S s 

U u 

MBAR M 

TIME t 


V v 

ASUM ^y i (MW) i 


PRESSURE p 

ENTHALPY h 


Distance along streamline from injection site 
divided by nose diameter 

Gas density, g/cm^ 

Droplet radius divided by initial droplet radius 

Debye length, cm 

Gas temperature, °K 

Compressibility factor 

Distance from injection site, cm 

Gas velocity, cm/sec 

Mean molecular weight of ions, g/mol 

Elapsed time after water injection, sec; 

<f\r 

t n = t n _ 1 + — ^x n - x n-1 j where n refers to values 
at end of forward integration step 

Droplet velocity, cm/sec 

Sum of mass fraction of species i, where 
14 

U5= I n 

i=10 

Gas pressure, dynes/cm^ 

Specific enthalpy of mixture, ergs/g 
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FORTRAN 

symbol 

Symbol 
in text 

Description 

WSUM 

Z w ‘ 

Sum of all gross reaction rates for all species i 


where i = 1 to 15, « 0, mol/cm^-sec 

EM 

M 

Mach number 

NSUBE 

N e 

Electron number density, electrons/cm^ 

FE 

Fe 

Efficiency factor for capture of electrons by 



water droplet 

F 

f 

Mass fraction of water evaporated 

CAPITALU 

U 

Normalized relative velocity between droplet and gas 

GAMMA 

n 

Specific concentration of species i after integration 



where i = 1 to 15, mol/g 

X 

Xi 

Mole fraction of species i where i = 1 to 15 

W 

Wi 

Net rate of production or disappearance of species i 



where i = 1 to 15, mol/cm^-sec 

LOWERCASE 

k i 

Reaction rate constant for reaction j where 

KSUBJ 

j = 1 to 45, cm/mol-sec 

UPPERCASE 

K i 

Equilibrium constant in mol/cm^ for reaction 

KSUBJ 

j = 1 to 6 and j = 17 to 21 and dimensionless for 
all other j (j = 1 to 40) 


AA,AB,AC 

k j( B ij " Gi j) 

Contribution of individual reactions to gross reac- 

etc., to BT 

tion W, mol-cmVg^-sec 


After each integration step, the concentrations are examined. If a concentration is 
negative, it is set equal to 0.0 and a comment "NEGATIVE GAMMAS" is printed. The 
concentrations are corrected to make the sum of the mass fractions of the species (ASUM) 
equal to 1.0. The corrected concentrations are printed when answers are printed. After 
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the addition of water and a second correction to make the sum of the mass fractions of 
the species equal to 1.0, the concentrations are again printed when answers are printed. 

Program Operation 

The following sketch shows the deck setup: 



The program stops are given in the following table: 


Stop 

Subroutine 

Reason 

Action 

11 

BASIC 

Iteration failed to converge 

Reexamine input data 

12 

CHECK 

Computing interval less than 10E-20 

l Vary ELE1 and ELE2 

13 

ERROR 

Input data in error 

[ Vary other parameters 
J (HEPS, Cl, CIMAX) 

100 

MAIN 

EOF, 5 

Normal stop 
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Flow Charts 

The flow charts for the D2331 main program and for the BASIC, CHECK, and ERROR 
subroutines are as follows: 


D2331 Main Program 
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Subroutine BASIC 


x2 

Used by CALINT to evaluate derivatives r , v, y i , i = 1 to 14 
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Subroutine ERROR 



Complete Program . 

The complete program including comments is reproduced in the following pages. 
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PROGRAM 023 21 ( INPUT , OUTPUT , TAP E5= I NPUT ,T APE6=0UTPUT ) 


STREAMTUBE // REPORT // PROGRAM 12-14-67 


STREAMTUBE ANALYSIS TO CALCULATE ELECTRON CONCENTRATION IN THE 
PRESENCE OF ACCELERATING ANO EVAPORATING WATER DROPLETS 

///// INPUT ///// 

IN = ONE CARD - 80 COLS - FOR IDENTIFICATION 

C$N AMI 

C RHOL = WATER DENS ITY ,G/CM**3 

C HL SPECIFIC ENTHALPY OF WATER, ERGS/G 

C DN REFERENCE LENGTH GIVEN AS NOSE DIAMETER, CM 

C EL LATENT HEAT OF WATER, ERGS/G 

C ENTHAL = INITIAL SPECIFIC ENTHA LP Y , ERGS /G 

C TW DROPLET SURFACE TEMPERATURE , DEG. K 

C FRG = DROPLET RECOVERY FACTOR FOR ZERO MASS TRANSFER 

C CPVW = HEAT CAPACITY OF VAPOR AT DROPLET SURF AC E , ER GS/ G-DEG. K 

C WASK = INITIAL MASS FLOW RATIO, WATER TO GAS 

C R INITIAL DROPLET RADIUS, CM 

C HASK = TOTAL SPECIFIC ENTHALPY OF SHOCKED GAS, ERGS/G 

C U2 INITIAL GAS VELOCITY AT INJECTION SITE, CM/SEC 

C P2 INITIAL GAS PRESSURE AT INJECTION SITE, DYNES/CM* *2 

C RH02 = INITIAL GAS DENSITY AT INJECTION S I TE , G/CM**3 

C T INITIAL TEMPERATURE OF MI XTURE , DEG.K 

C PRDEL = PRINT ANSWERS EVERY PRDEL X. =0.1 UNLESS INPUT OTHERWISE 

C 

C VAR ( 1 ) = X 

C VAR ( 2 ) = R HAT SQUA RED , ( DROPL ET RAD IUS/ IN ITI AL DROPLET RADIUS>**2 

C VAR ( 3 ) = V, DROPLET VELOC I TY , CM/S EC 

C V AR ( 4 ) = GAMMA FOR 1=1 SPECIFIC CONCENTRATION OF H20,M0LE/G 

C VAR ( 5 ) = 1=2 02 

C V AR ( 6 ) = 1 = 3 N2 

C 

C ETC FOR 14 GAMMAS ETC 

C 

C 0ERI1) = DERIVATIVE OF R HAT SQUARED 

C DER ( 2 )= V DOT 

C DER ( 3 ) = GAMMA ( 1 ) DOT 


100000 

200000 

300000 

400000 

500000 

600000 

700000 

800000 

900000 

1000000 

1100000 

1200000 

1300000 

1400000 

1500000 

1600000 

1700000 

1800000 

1900000 

2000000 

2100000 

2200000 

2300000 

2400000 

2500000 

2600000 

2700000 

2800000 

2900000 

3000000 

31C0000 

3200000 

3300000 

3400000 

3500000 

3600000 

3700000 

3800000 

3900000 

4000000 

4100000 

4200000 


09 

CD 
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o 


000003 

000003 

000003 

000003 

000003 

000003 

C00003 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c$ 

c 

c 

c 

c 

c 

c 

c 


DER{4)= G A MMA ( 2 ) DOT 
ETC FOR 14 GAMMAS 

ELE 1 = AN ARRAY OF 16 VALUES USED BY THE INTEGRATION SCHEME TO 

CONTROL THE SIZE OF THE COMPUTING INTERVAL 
0. 0.LT.ELE1.LE.65.0 

INITIALLY 0.5 UNLESS INPUT OTHERWISE 
ELE2 = AN ARRAY OF 16 VALUES USED BY THE INTEGRATION SCHEME TO 
CONTROL THE SIZE OF THE COMPUTING INTERVAL. 

E LE2.LT.ELE1 

INITIALLY 0.1 UNLESS INPUT OTHERWISE 
Cl = INITIAL COMPUTING I NT ERV AL , = . 00 01 2 20 703 1 2 5 UNLESS INPUT 

CIMAX = MAXIMUM COMPUTING INTERVAL , = .1 UNLESS INPUT OTHERWISE 
XTERM = TERMINAL VALUE OF X 

HEPS = PERCENT BY WHICH ENTHALPY ALLOWED TO VARY IN TWO 

SUCCESSIVE STEPS WHEN ENTHA LPY . GT . 50000 . .OR . ENTHALPY 
•LT. -500000. TEST IN SUBROUTINE CHECK 

0.1 UNLESS INPUT OTHERWISE 

FE = EFFICIENCY FACTOR FOR CAPTURE OF ELECTRONS BY WATER 

DROPLET OR CATALYTIC EFFICIENCY 

NTLUP = NUMBER OF POINTS IN VARIX AND VARDP TABLES 
VARIX = TABLE OF X VALUES FOR P/P2 

VARDP = TABLE OF P/P2 VALUES CORRESPONDING TO VARIX 
I BUG = 0 NO DEBUG PRINT OUT, 0 UNLESS INPUT 

= 1 FOR DEBUG PRINT OUT AT DESIRED PLACES 
ARRAY OF 20 


STOPS 

STOP 11 IN BASIC AT ITR1 
STOP 12 IN CHECK WHEN C I .LT . 1 . OE-2 0 
STOP 13 IN ERROR WHEN SOME INPUT IN ERROR 
STOP 100 EOF, 5 NORMAL STOP 

COMMON VAR (20) , DER ( 2 0 ) , ELE 1 ( 20 ) , ELE2 ( 20 ) , Cl , SP EC , NDE , ELT , NE RR 
COMMON RHQL , HL , ON, EL , TW ,CDC ,COFM , FRO , CP VW, WASK , R , HASK 

COMMON T 

COMMON' El, E2,P2,U2,RH02,V2 

COMMON F,NE» Z,XI( 15) , FAC 1 , F AC2 , U , RHOU ,RH02 U2 , R HO 
COMMON KLC145) ,KUC(40) , MB AR , FE , C APU, L AMD, P 

COMMON AA, AB, AC, AD, AE , AF , AG, AH , A I , A J , AKK , AL , AM , AN , AO, AP , AQ , AR , 

1 AS,AT,AU,AV, AW, AX, AY, AZ,BA,BB,BC, BD , BE , BF , BG , BH, B I , B J , BKK, 

2 BL,BM,BN,BO,BP,BQ,BR,BS,BT 


4300000 

4400000 

4500000 

4600000 

4700000 

4800000 

4900000 

5000000 

5100000 

5200000 

5300000 

5400000 

5500000 

5600000 

5700000 

5800000 

5900000 

6000000 

6100000 

6200000 

6300000 

6400000 

6500000 

6600000 

6700000 

6800000 

6900000 

7000000 

7100000 

7200000 

7300000 

7400000 

7400001 

7500000 

7600000 

7700000 

7800000 

7900000 

8000000 

8100000 

8200000 

8300000 

8400000 
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000003 

000003 

000003 


000003 


000003 


000003 


000003 

000003 

000003 

000003 

000003 

000003 

000003 

000003 


000003 


c 

c 

c 

c 

c 

c 


c 


c 

c 

c 


c 

000003 

000003 


000003 


COMMON NTL UP » VARI X ( 40 ) , VARDPI 40 ) 

EQUIVALENCE ( A A, ALPDUM ( 1 ) ) 

DIMENSION ALPDUMI46) 

LABI IN MAIN, BASIC, FOFX CHECK 
COMMON/LABI /I BUG, CU VAR, A1 ,A2, A3, A4,A5 , A6 , TL I M 1 , TL I M2 , TL I M3 
LAB 2 IN MAIN, BASIC, CHECK, FOFX 
COMMON/ LAB 2 /ENTHAL 

LAB3 IN MAIN , BASIC AND CHECK 
COMMON /LAB3/ TPRE V ,KEY, HEP S , KA Y 
LAB4 IN MAIN AND BASIC 

COMMON /LAB 4/ MW I 15 ) , W ( 1 5 ) , GAMM A , HC ITR , B , DELTX , EM 

DIMENSION All 15,4) ,A2( 15,4) t A3 ( 15,4), A4( 15,4) ,A5( 15,4) ,A6( 15, 4) 
DIMENSION CUVARI20) , IBUGI20) 

DIMENSION SPECIE (15), FGAM( 15 ) 

DIMENSION IN ( 8 ) 

REAL I N I TX » MVI 
REAL LH, NE 

DOUBLE PREC ISION VAR 

MOLECULAR WEIGHT OF SPECIES 

DATA (MW (I), 1=1, 15) / 18. ,32. ,28. ,2. ,30. ,17. , 14. , 16. ,1 .,30. ,28., 
1 14. , 1 6. ,1 . , 5.486E-4/ 

THE 15 SPECIES 

DATA { SP EC I E ( I ) ,1 = 1,15) / 3HH20 , 2H02 , 2HN2 , 2HH2 , 2HN0 , 2H0H , 1HN, 1H0, 
1 1HH , 3HN0+ , 3HN2+,2HN+,2H0+,2HH+,2HE-/ 

DATA (Al(I) ,1=1,60)/ 4.20307985,3.81416055,3.69148176,3.45380035, 

1 3.87682186,3.44021610, 2.5 ,3.04980070,2.50020545,3.68482747, 

2 3.76510471,2.80887381,2.50047586, 2.5, 2.5, 

3 3.48589794,3. 85044317,3. 34031165,2.97106943,3.71537602, 

4 3.01845824,2. 75947630,2.58926700,2.94053868,3.26815342, 

5 3.09030675,2.55758547,2.54621952,2.5 ,2.5 , 30*0./ 

DATA (A2(I ) ,1 = 1,60) / -2.14625506E-3,-3. 2409 19 83E-3, 

1 -1.33248199E-3, 4 . 2 1 2 22704E-4 , -3 . 147 87008E-3 , 5. 6 7C61 0 81 E-4, 0., 

2 — 2. 3867 48 35E-3,— 1.36545791 E- 6,-1. 27472750E-3,— 1.98295059 E- 3, 

3 -1.48986257E-3,-2.872796ClE-6, 0.,0., 

4 1.8 17201 2 IE-3, 3 . 899 89763E-4 , 7 . 7193605 IE-4 , 7. 60556499E-4, 

5 4.81 4032 7 8E-4 , 7. 9 16 17340E-4.-3 . 58661 047E-4 ,- 1 . 1 1816431 E-4, 

6 -6. 7 006204 4E -4, 8 . 74276324E-4 , 1 *03070 1 5 8E-3 ,-7. 929528 60E-5, 

7 -5.12662872E-5, 0., 0., 30*0./ 


8500000 

8600000 

8700000 

8800000 

8900000 

9000000 

9100000 

9200000 

9300000 

9400000 

9500000 

9600000 

9700000 

9800000 

9900000 

10000000 

10100000 

10200000 

10300000 

10400000 

10500000 

10600000 

10700000 

10800000 

10900000 

11000000 

11100000 

11200000 

11300000 

11400000 

11500000 

11600000 

11700000 

11800000 

11900000 

12000000 

12100000 

12200000 

12300000 

12400000 

12500000 

12600000 

12700000 
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to 


000003 


000003 


000003 


000003 


000003 


c 

c 


c 

c 

c 

c 

C00003 


DATA ( A3 ( I ) ,1 = 1,60) / 6. 72040236E-6, 1 . 04955 3 52 E-5 , 

1 2.64950496E-6.-1. 10268215E— 6, 8 . 40430066E-6 , -1 . 79909482E-6 , 0., 

2 4. 3187221 7E-6, 3 . 0 1 1737 18 E-9 , 2 . 4846 1737 E-6, 4. 53722 826E-6, 

3 2. 8757776 0E-6, 5 . 66990999E-9 , 0.» 0., 

4 -3.41338632E-7,-2.89418093E-8,-l . 7562 2629E-7 , - 1 . 16571433E-7, 

5 -9. 721942 39E- 8,-1 . 390 467 62 E— 7 , 1 . 23? 19036E-7 , 4. 3605 8822E-8 , 

6 2.93518521E-7.-2. 1 84678 1 8E-7 1 . 559 1 0644E-7 , 3 . 1 1 1 52083E- 8, 

7 6.71501704E-9, 0., 0., 30*0./ 

DATA ( A4( I) ,1=1,60) /- 5. 56262948E-09 1 . 03969441 E-08 , 

1 -9.7842 06 15E-10, 1 . 2 2466 55 5E-9 , -7 . 20675049E-C 9 , 2 . 3209 1 070E-09 , 0 , 

2 -3.615676 0 2E-09,-2.745164 50E-12,-8.1 1 027 164E-1 0,-2 . 9 130627 IE-09, 

3 -2. 5121 89 7 4E- 09, -5 . 16787546E-12, 0., 0., 

4 2.68468828E-11, 1 . 1607 683 9E-1 2 , 1 . 6 14 16 590E- 1 1 , 8 .39248279E-1 2 , 

5 9.43 59166 IE- 12 , 1 . 1 2 28 27 53E-1 1 , - 1 . 08890 1 56E- 1 1 , -4 . 77729929E-1 2 , 

6 -4.67055452E-11, 2 . 271 1 6436E-11 , 9 .92335 507E-1 2, - 2. 9 281633 5E-1 2, 

7 1 . 653499 3 7E- 12 , 0., 0., 30*0./ 

DATA (A5(I ) ,1 = 1,60) / 1 . 7433 893 8E-12 , 3 . 5 22 80537E- 12 , 

1 -9. 817897 75E - 14, -3. 64412032 E- 13, 2 . 1 536261 2E- 1 2 , - 8. 4944147 IE-1 3 , 
20,1. 1496373 0E- 12, 8 .9496 199 7E-16 , - 1 . 55 190996E- 13 , 5 .75631866E-1 3 , 

3 8. 2272882 7E -13, 1 . 87850329E- 1 5 , 0., 0., 

4 -7. 4634477 6E- 16 , -3 . 2309998 2 E-17 , -4. 45648938E- 16, - 2 . 1 2 1 01689E- 1 6, 

5 -3 . 0279653 3E-16, -3. 1798705 7E-16, 3. 06742425E-16, 1.68389768E-16, 

6 2.41598190E-15,-7.03854576E-16,-2.28531083E-16, 8 .6969 1 886E-1 7, 

7 -1. 0421302 4E -16, 0., 0., 30* 0./ 

DATA (A6(I) ,1 = 1,60) / -2 . 8749849 1 E4,-2 . 193 1860 5E1 , -1 . 907136 53E 1 , 

1 4.18216914 , 1.07744359E4,4.68483226E3,5.66215000E4, 

2 2 . 9649358 0E4 , 2 . 59754478E4, 1 . 1 8655 57 IE 5 , 1 . 807793 86E 5 , 

3 2 . 2 5 28108 7E 5 , 1 . 87666445E 5 , 1 . 83 7550 00E5 , 0., 

4 -2.87104401E4,-2.61240243E2,-4.75889275E1, 2 . 201 8 8265E 2 , 

5 1 . 0580474 2E4 , 4. 8 48 50446E 3 , 5 . 64949247E4, 2 . 9724623 1 E4 , 

6 2. 5782397 8E4, 1 . 1 865893 5E 5 , 1 . 8 0896602E 5 , 2 . 253 1 2369 E5 , 

7 1. 8764322 4E 5, 1 . 8 375 50 00E 5 , 0., 30*0./ 

INPUT 

NAME LIST / NAM1/ RHOL , HL , ON , EL , ENT HAL , TW , FRO , C PVW , 

1 MASK »R,HASK* IBUG, U2, P2,RH02,T , 

2 PRDEL , VAR,DER,ELE1 ,ELE2»CI ,C1MAX,XTERM, HEPS , 

3 FE,NTLUP,VARI X, VARDP 

***** SET NUMERICAL CONSTANTS ***** 

TL IM1 = 1000 . DEG.K 
TL I Ml=l 000 . 


12800000 

12900000 

13000000 

13100000 

13200000 

13300000 

13400000 

13500000 

13600000 

13700000 

13800000 

13900000 

14000000 

14100000 

14200000 

14300000 

14400000 

14500000 

14600000 

14700000 

14800000 

14900000 

15C00000 

15100000 

15200000 

15300000 

15400000 

15500000 

15600000 

15700000 

15800000 

15900000 

16000000 

16100000 

16200000 

16300000 

16400000 

16500000 

16600000 

16700000 

16800000 

16900000 

17000000 
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c 


TL I M2= 1 5000 . DEG K 

17100000 

C00005 



TLX M2 = l 5000 . 

17200000 


c 


TLI M3=1000000. 

17300000 

000006 



TLI M3= 10000 00 . 

17400000 

000010 



SAVECI=. 000 12207031 25 

17500000 

000011 



PRDEL = . 1 

17600000 

000013 



C I M AX= . 1 

17700000 

000014 



FAC 1= 4./ ( 3. *2.7182818284 * 3. 141 5926536**. 5 ) 

17800000 

000022 



FAC2= l./I 2. *2. 7182818284**1. 5) 

17900000 


c 


*** READ IDENTIFICATION CARD 

18000000 

000027 



RE AD ( 5 * 2 J ( INI I), 1=1, 8) 

18100000 

000040 


2 

FORMAT! 8A1 0) 

18200000 


c 


KEY SET = 1 IN CHECK TO ALLOW TESTING OF ENTHALPY AFTER FIRST INTERVAL 18300000 


c 


KAY SET = 1 IN BASIC WHEN T.LT.O. 

18400000 


c 


BEGIN NEW CASE AT EFN 1 

18500000 


c 


*** INITIALIZE 

18600000 

000040 


1 

CONTINUE 

18700000 

000040 



HEP S = 0 . 1 

18800000 

000042 



KEY=0 

18900000 

000043 



C I = SA VEC I 

19000000 

000044 



KAY = 0 

19100000 

000045 



DO 4 1=1,46 

19200000 

000046 



ALPDUM ( I 1=0 . 

19300000 

000047 



I F I I . GT . 20 ) GO TO 4 

19400000 

000053 



VARI I 1=0. 

19500000 

000055 



DERI I > = 0. 

19600000 

000056 



CUV AR l I ) =0 

19700000 

000057 



I BUG! I 1=0. 

19800000 

000060 



ELE1I I)= .5 

19900000 

000061 



ELE2I I )= .1 

20000000 

000063 



I F I I . GT .15) GO TO 4 

20100000 

000066 



W( 11=0. 

20200000 

000067 


4 

CONTINUE 

20300000 


c 



20400000 


c 


*** READ NAMELIST INPUT 

20500000 

000071 



READ! 5 , NAM1 ) 

20600000 


c 



20600001 

000074 



I F I EOF , 5 ) 7,9 

20600002 

000077 


7 

PRINT 8 

20600003 

000103 


8 

FORMAT!// 13H EOF, 5 - NAMELIST//) 

20600004 

000103 



STOP 100 

20600005 

000105 


9 

CONTINUE 

20600006 


c 



20600007 


4 * 

CO 
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000105 


WRITE(6,N4M1) 

20700000 


C 


0012 20800000 


C 

***** COMPUTE CASE CONSTANTS FROM INPUT 

20900000 


c 


21000000 

000110 


I N I TX= V AR( 1 ) 

21100000 

000112 


T I ME = 0 . 

21200000 


c 

EQUATION ( 10 ) 

21300000 

000113 


VAR I 3 ) = I ( 600. *RH02*U2*WASK)/3. 141 59 26535 89793 ) / 

21400000 



1 I RHOL*VAR ( 2) *SQRT ( VAR ( 2 ) ) ) +1. 

21500000 

000151 


V2=V AR ( 3 ) 

21600000 

000153 


U=U2 

21700000 

000154 


RH02U2=RH02 *U2 

21800000 

000155 


RH0U= RH02U2 

21900000 


c 

*** PRINT INITIAL VALUES *** 

22000000 

000156 


PRINT 5 

22100000 

000162 


5 FORMAT ( 1X,*D2331 - STRE AMTUBE / /REPORT / / PROGRAM*/) 

22200000 

000162 


PRINT 3, ( INI I ) ,1 = 1,8) 

22300000 

000174 


3 FORMAT! 8A 1 0/ ) 

22400000 

000174 


PRINT 6, WASK,R,0N,P2,RH02,ENTHAL,T,U2,V2,FE 

22500000 

000224 


6 FORMAT ( 26H INJECTION PARAMETERS-—/ 15X 7HW STAR=E1 5 . 8, 6X 

22600000 



1 2HR= El 5.8, 10X 3HDN=E15.3/ 

22700000 



2 22H INITIAL VALUES — P=E15.8,8H RHO=E15.8, 

22800000 



213H ENTHALPY=E15.8/22H T=E15.8,8H U2= 

22900000 



3E 15 . 8 , 1 3H V2=E 15.8/ 19X 3HFE=E15.8) 

23000000 

000224 


DO 10 1=1,15 

23100000 

00 C 226 


10 PRINT 11 ,SPECI E( I ) , I 

23200000 

000237 


11 FORMAT (2X, A10, 15 ) 

23300000 

000237 


F=0. 

23400000 

000240 


XPREV=VAR( 1 ) 

23500000 

000242 


TPREV=T 

23600000 

000244 


DO 42 1=1,14 

23700000 

000245 


42 FGAMI I ) = VAR ( l +3 ) 

23800000 

000252 


DIG I 2= 16.*(FGAM(1)+2.*FGAM(2)+FGAM{5)+FGAM(6)+FGAM(8)+FGAM(10) 

23900000 



1 +FGAMI 13) ) 

24000000 

000263 


D2G 12= 14.* (2. *FGAM( 3 ) +FGAM l 5 ) + FGAMI 7 ) +FGAMI 1 0 ) +2 . *FG AM { 1 1 ) 

24100000 



1 +FGAMI12) I 

24200000 

000273 


DO 40 1=1,17 

24300000 

000274 


40 CUVAR ( I ) =V AR I I ) 

24400000 

000301 


NDE=1 6 

24500000 

000302 


PHMAX=65 . 

24600000 

000303 


XPRINT=VAR (1) 

24700000 

000305 


I F I N = 0 

24800000 

000306 


IW AY= 1 

24900000 


APPENDIX 



000307 

000313 


000313 

000313 

000314 

000315 


000315 

000331 

000334 

000336 


C *** DEFINE OUTPUT 

PRINT 45 

45 FORMAT (/5X, 61 HX=D I ST ANC E ALONG STREAMLINE FROM INJECTION SITE/NOSE 

1 DI A ME TER/ 5 X, 23HRH0 = GAS DENSITY , G/CM**3/ 5X , 54HR HAT^DROPLET RAOIUS 

2 DIVIDED BY INITIAL DROPLET R AD I US / 5X , 1 09HL AMBD A SUB D=DEBYE LENGT 
3H, CHARACTER 1ST IC DISTANCE OF FIELD BEYOND WHICH THE EFFECT OF A CH 
4ARGE IS NOT FELT, CM// 

45 X, 33HTEMPERATURE=GAS TEMPERATURE ,DEG.K, 27X,21HU=GAS VELOCI TY , CM/S 
5EC/5X , 2 4HZ= COMPRESS I B I LI TY F ACTOR ,36X , 43HM BAR = MEAN MOLECULAR WEIG 
6HT FOR IONS ,G/MOLE/ 5X» 50HS=DI STANCE ALONG STREAMLINE FROM INJECTIO 
7N SITE, CM, 10X,43HTI ME=EL APSED TIME AFTER WATER INJECT ION , SEC// 
85X,25HV=DR0PLET VELOC I T Y, CM/SEC , 35X , 39 HENTHALP Y=SPECI F I C ENTHALPY 
9 OF GAS , ERG/ G/ 5X ,38HA SUM=SUM OF MASS FRACTIONS OF SPECIES, 22X, 
A61HW SUM=SUM OF ALL THE GROSS REACTION RATES, MOLE/ CM** 3 -SEC ,=0.0/ 
B5X,33HPRESSURE=GAS PRESS UR E, DYNES/CM* *2, 27X, 1 1HMACH NUMBER// 

C 5X , 46HN SUB E = ELECT RON CONCENTRATION, ELECTRONS /CM* *3/ 

D5X, 86HFE=EFFI CIENCY FACTOR FOR CAPTURE OF ELECTRONS BY WATER DROPL 
EET OR CATALYTIC EFFICIENCY// 

F5X, 30HF=FRACT I ON OF WATER EVAPORATED/ 

G5X , 62HC API T AL U=NORMAL I ZED RELATIVE VELOCITY BETWEEN DROPLET AND G 
HAS/ ) 

C 
C 
C 

101 CONTINUE 
ELT= 0 
SPEC=0 

102 CONTINUE 
C 
c 
c 
c 
c 
c 
c 

CALL CALINT ( V AR , DER , EL El , EL E2 , C I , SP EC ,N DE , CUVAR , E L T ,C I M AX , NERR , 

1 PHMAX ) 

IWAY=-1 *IWA Y 

C GO TO 51 TO PRINT INITIA3 CONDITIONS 

IF{ INITX.EQ.VAR(l)) GO TO 51 
IF( IWAY.EQ.-l. AND . V ARI 1) .NE.INITX) GO TO 102 
C 

C 5/13/67 TREAT THE PRODUCT RHO*U AS CONSTANT OVER AN INTERVAL 


REINITIALIZE EACH COMPUTING INTERVAL SINCE GAMMAS ARE 
ALTERED 


CALINT EVALUATES DERIVATIVES IN SUBROUTINE BASIC 
SEE SUBROUTINE BASIC FOR DETAILS OF CALCULATION PROCEDURE 
FOR R HAT SQUARED DOT, V DOT AND GAMMA 1 DOT THRU GAMMA 
14 DOT AND THUS THE EVALUATION OF THE BULK OF THE 
EQUATIONS AND AUXILIARY RELATIONSHIPS 
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25100000 

25200000 

25300000 

25400000 

25500000 

25600000 

25700000 

25800000 

25900000 

26000000 

26100000 

26200000 

26300000 

26400000 

26500000 

26600000 

26700000 

26800000 

26900000 

27000000 

27100000 

27200000 

27300000 

27400000 

27500C00 

27600000 

27700000 

27800000 

27900000 

28000000 

28100000 

28200000 

28300000 

28400000 

28500000 

28600000 

28700000 

28800000 

28900000 

29000000 

29100000 

29200000 


oi 
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03 


000346 



IF ( WASK. EQ . 0. ) VAR ( 3 ) =U 

000351 



RH0U= RH0*U 

000353 

r 

51 

CONTINUE 

000353 

u 


00 59 1=1,14 

000355 


59 

FGAM ( I )=VAR ( I +3 » 

000362 


57 

CONTINUE 


c 


*** COMPUTE ASUM,WSUM,GAMMA<15) ,X,S,Z,TIME 


c 


ASUM SHOULD BE = 1.0 

000362 



ASUM= 18.4 FGAM(1)+32.*FGAM(2)+28.*<FGAM<3)+FGAM(11)) + 2.*FGAM(4) 



1 

+30.*(FGAM (5 ) +FGAM (10)1 + 17.*FGAM(6> + 14 . 4 ( FGAM ( 7) +FGAMI 12 )) 


r 

2 

! +16. *( FGAM (8 ) +FGAM ( 13 )) + FGAM ( 9 ) + FGAMI14) 


c 


WSUM SHOULD BE = 0.0 1=1,14 

000414 


63 

CONTINUE 

000414 



WSUM=0 . 

000415 



DO 65 1=1,14 

000417 


65 

HSUM= WSUM + MW( I )* W ( I ) 

000424 



FWASK= F*WASK 

000426 



FGAMI 15)=0. 

000427 



DO 68 1= 10,14 

000430 


68 

FGAMI 151= FGAMI 15)+FGAM( I ) 

000434 



GAMMA=0. 

000435 



DO 67 1=1,15 

000436 


67 

GAMMA=GAMM A+FGAMI I ) 

000442 



DO 66 1=1,15 


c 


EQUATION (39) 

000443 


66 

X I l I ) = FGAM I I ) /GAMMA 

000447 



S= VAR I 1 ) *D N 


c 


EQUATION 141) 

000455 



Z = GAMMA*! 1.0+ FWASKJ/I .034674063 + FWASK/18.) 

000462 



IF I VARI3 ) • EQ. 0 . ) TIME =0. 

000464 



IF I VAR I 3) . EQ.O . ) GO TO 69 

000465 



TIME=TIME+DN*( VARI 1 l-XPREVJ /VAR I 3 ) 

000513 


69 

XPRE V= VARI 1 ) 


c 


TPREV USED IN CHECK 

000515 



TPRE V=T 


c 


*** TERMINATION TEST 

000517 



IFIVARID.LT. XTERM) GO TO 52 

000521 



IF! IF IN.NE . 0) GO TO 82 

000522 



XPRINT=VAR I 1) 

000524 



I F I N= 1 

000525 



GO TO 54 
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29400000 

29500000 

29600000 

29700000 

29800000 

29900000 

30000000 

30100000 

30200000 

30300000 

30400000 

30500000 

30600000 

30700000 

30800000 

30900000 

31000000 

31100000 

31200000 

31300000 

31400000 

31500000 

31600000 

31700000 

31800000 

31900000 

32000000 

32100000 

32200000 

32300000 

32400000 

32500000 

32600000 

32700000 

32800000 

32900000 

33000000 

33100000 

33200000 

33300000 

33400000 

33500000 


> 

d 

d 

fd 


2 

d 





000526 


52 CONTINUE 

33600000 


C 

*** TEST PRINT FREQUENCY 

33700000 

000526 


IF! XPRINT-VAR! 1 ) .GT.0.0) GO TO 81 

33800000 


C 

PRINT ANSWERS EVERY PRDEL IN X 

33900000 

000532 


XPR I NT= XPR I NT+PRDEL 

34000000 

000533 


54 CONTINUE 

34100000 

000533 


I KY=0 

34200000 


C 

*** PRINT ANSWERS 

34300000 

000534 


PRINT 70 

34400000 

000540 


70 FORMAT! 5X 1HX 10X 3 HR HO 14X 11HTEMPERATURE 6X 1HU 16X 1HV 16X 

34500000 



1 8HENTHALPY 9X 7HN SUB E 10X 1HF/ 

34600000 



1 16X 5HR HAT 12X 1HZ 16X 5HM BAR 12X 5HA SUM 12X 5HW SUM 12X 2HFE 

34700000 



215X 9HCAPIT AL U/ 15X 1 2HLAMBDA SUB D 6X 1HS 16X 4HTIME 

34800000 



3 13X 8H PRES SURE 9X 2HEM/) 

34900000 


C 


35000000 

000540 


F X = VAR ( 1 ) 

35100000 

000542 


FV=VAR! 3 ) 

35200000 

000544 


FRHAT= SQRT ( V AR ( 2 } ) 

35300000 

000547 


PRINT 71 » F X ,RHO,T,U,FV ,ENTHAL,NE,F 

35400000 

000572 


71 FORMAT (Fll .6.7E17.8 ) 

35500000 


c 


35600000 

000572 


PRINT 72 , FRHAT , Z , MB AR , ASUM , WSUM ,FE ,C AP U 

35700000 

000614 


72 FORMAT1E28. 8.6E17.8) 

35800000 


c 


35900000 

000614 


PRINT 73, LAMO,S,TIME,P,EM 

36000000 

000632 


73 F0RMAT!E28. 8,4E17.8/) 

36100000 


c 


36200000 

000632 


PRINT 74, ( FGAM! I) , 1 = 1,15) 

36300000 

000644 


74 F0RMAT(11X,74HGAMMA, SPECIFIC CONCENTRATION OF SPEC I ES , MOLES /G ,0N R 

36400000 



1ETURN FROM INT EGRAT I ON/ 1 2X , 1 HI , E 1 5 . 8 , 4E 1 7 . 8/ 1 2X , 1H 6 ,E 1 5 . 8 ,4E1 7 . 8/ 

36500000 



211X,2HU,E15.8,4E17.8/) 

36600000 


c 


36700000 

000644 


PRINT 75, (XI ( I) ,1=1,15) 

36800000 

000656 


75 FORMAT ( 11X 25HX.MOLE FRACTION OF SPECIE/ 3 ( E28 . 8 , 4E 17 . 8/ ) ) 

36900000 


c 


37000000 

000656 


PRINT 76, ( Will ,1=1 ,15 ) 

37100000 

000670 


76 FORMAT! 11X, 67HW.NET RATE OF PRODUCTION OR DISAPPEARANCE OF SPECIES 

37200000 



1,M0LE/CM**3-SEC/12X,1H1,E15 .8 , 4E1 7. 8/ 12 X, 1H6 , E 15. 8 ,4E 17 . 8/ 1 IX, 2H1 1 

37300000 



2.E15.8, 4E17.8 /) 

37400000 


c 


37500000 

000670 


PRINT 77, ! KLC ( I) ,1=1,45) 

37600000 

000702 


77 FORMAT! 11X.90HLOWER CASE K SUB J , J = 1 , 45 , REACT I ON RATE CONSTANTS FO 

37700000 



1R A PARTICULAR REACTION ,CM**3/M0LE-SEC/6 ( E28.8.6E17.8/ JE28.8, 2E17. 

37800000 


-a 


APPENDIX 




28/) 


37900000 


C 



38000000 

000702 


PRINT 78, ( KUC ( 1 ) , 1=1 ,40 ) 


38100000 

000714 


78 FORMAT! 11X ,76HUPPER CASE K SUB J , J=1 , 40 , EQU I L I BR I UM CON STANTS , MOLE 

38200000 



1S/CM**3 OR DIMENSIONLESS/ 5 ( E2 8 . 8 , 6E 17 . 8/ ) , E 28 . 8 , 4E 17. 8/ ) 


38300000 


CCCCC 


38400000 

000714 


PRINT 179 


38500000 

000720 


179 FORMAT! 11X,76HC0NTR IBUTION CF INDIVIDUAL REACTIONS TO GROSS REACTI 

38600000 



ION, W, MOLE- CM**3/G** 2-SEC) 


38700000 

000720 


PRINT 79 


38800000 

000724 


79 FORMAT ( 16X 2HA A 15X 2HAB 15X 2HAC 15X 2HAD 1 5X 2HAE 1 5X2H AF 

15X2HAG/ 

38900000 



1 16X 2HAH 1 5 X 2HAI 15X 2HAJ 15X 2HAK 15X 2HAL 15X 2KAM 15X 

2HAN/ 

39000000 



2 16X 2HA0 1 5 X 2 HAP 15X 2HAQ 15X 2 HAR 15X 2 HAS 15X 2HAT 15X 

2HAU/ 

39100000 



3 16X 2HAV 1 5 X 2 HAW 15X 2HAX 15X 2HAY 15X 2HAZ 15X 2HBA 15X 

2HBB/ 

39200000 



4 16X 2HBC 15 X 2HBD 15X 2HBE 15X 2HBF 15X 2HBG 15X 2HBH 15X 

2HBI / 

39300000 



5 16X 2HBJ 1 5 X 2HBK 15X 2HBL 15X 2HBM 15X 2HBN 15X 2HB0 15X 

2HBP/ 

39400000 



616X 2HBQ 1 5 X 2HBR 1 5X 2HBS 15X 2HBT/ ) 


39500000 


C 



39600000 

000724 


80 F0RMAT!6!E28.8 ,6E17 .8/ ) ,E28.8,3E17.8/) 


39700000 

000724 


PRINT 80, AA, AB, AC, AD, AE,AF,AG, AH, A I , A J , AKK , AL , AM , AN , 


39800000 



1 AO, AP , AQ , AR , AS, AT , AU , AV , AW , AX , AY , A Z , B A , BB , 


39900000 



2 BC,BD,BE,BF,BG,BH,BI, B J , BKK , BL , BM , BN , BO , BP , 


40000000 



3 BQ,BR,BS,BT 


40100000 


C 



40200000 

001064 


81 CONTINUE 


40300000 


C 

*** TEST FOR NEGATIVE GAMMAS, IF NEGATIVE SET = 0.0 


40400000 

001064 


DO 61 1=1,15 


40500000 

001066 


IFtFGAM! I) . GE.O.O) GO TO 61 


40600000 

001070 


PRINT 37 


40700000 

001073 


37 FORMAT ! 19H NEGATIVE GAMMAS) 


40800000 

001073 


V AR ( I +3 )=0 . 


40900000 

001076 


FGAM ( I ) =0. 


41000000 

001077 


61 CONTINUE 


41100000 


c 



41200000 


c 

REPLACE GAMMAS WITH ADJUSTED VALUES 


41300000 


c 



4140000C 

001101 


I P ATH= 1 


41500000 

001102 


30 CONTINUE 


41600000 

001102 


B 1= ( D 1G I 2 +(16./18.>* FWASK)/(1.0+FWASK) 


41700000 

001111 


B2= D2GI 2/ ( 1.0+FWASK) 


41800000 

001113 


83= ( <2./18.)*FWASK)/(1.0+FWASK) 


41900000 

001120 


Dl= 81/(16.*! FGAM! 1 ) +2 . *FG AM ( 2 ) +FG AM ( 5 ) +FGAM ( 6 ) +F GAM < 8 ) +FG AM ( 1 0 ) 

42000000 



1 +FGAMU3) ) ) 


42100000 
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001132 


02= B2/(14.*(2.*FGAM(3J+FGAM(5)+FGAM<7)+FGAM< 10 ) +2 .*FGAM( 11) 

42200000 



1 +FGAM( 12) ) ) 

42300000 

001143 


D3= B3/(2.*FGAM(1 )+2.*FGAM(4)+FGAM(6) +FGAM( 91+FGAM (14) ) 

42400000 

001152 


IF( FGAM( 1) . EQ. 0. .AND. FGAM ( 4 ) . EQ . 0. .AND. 

42500000 



1 FGAMI6) . EQ.O. .AND. FG AMI 9 ) . EQ .0 . .AND. F GAM ( 14) . EQ. 0. ) D3=0. 

42600000 

001167 


D4= (2.*D3+16.*D1)/18. 

42700000 

001174 


05= (14. *02+16. *Dl)/30. 

42800000 

001177 


D6= ( D3+16. *D1 > / 17 . 

42900000 


C 


43000000 


C 

EQUATION (32) 

43100000 

001202 


FGAM( 2 )= FG AM ( 2) *D1 

43200000 

001203 


FG AM ( 8 ) = FGAM( 8 ) *D1 

43300000 

001204 


FGAM( 13)=FGAM( 13)*D1 

43400000 


C 


43500000 


c 

EQUATION (33) 

43600000 

001205 


FGAM ( 3 ) = FGAM( 3) * D2 

43700000 

001207 


FG AM ( 7 ) = FGAM( 7) * D2 

43800000 

001210 


FGAM4 11)=FGAM( ID* D2 

43900000 

001211 


FGAM( 12 )=FG AM ( 12 )* D2 

44000000 


c 


44100000 


c 

EQUATION (34) 

44200000 

001212 


FGAM ( 4 ) =FG AM ( 4 ) *D3 

44300000 

001213 


FGAM( 9 )=F G AM ( 9 ) *D3 

44400000 

001214 


FGAM( 14)=FGAM( 14) *D3 

44500000 


c 


44600000 


c 

EQUATION (35) 

44700000 

001215 


FGAM( 1 )=FGAM( 1 )*D4 

44800000 


c 


44900000 


c 

EQUATION (36) 

45000000 

001217 


FGAM( 5 )=FG AM ( 5 )* D5 

45100000 

001220 


FGAM ( 10 ) =FG AM ( 10)* D5 

45200000 


c 


45300000 


c 

EQUATION (37) 

45400000 

001221 


FGAM ( 6 ) = FGAM( 6) * D6 

45500000 


c 

ASUM SHOULD BE = 1.0 

45600000 

001223 


ASUM= 18.* FGAM(1)+32.*FGAM(2)+28.*(FGAM(3)+FGAM( ID) + 2.*FGAM(4) 

45700000 



1 +30.*(FGAM(5)+FGAM(10) ) + 17.*FGAM(6) + 14. * ( FGAM ( 7) +FGAM( 1 2 ) ) 

45800000 



2 +16.*(FGAM(8)+FGAM(13) ) + FGAM ( 9 ) + FGAM( 14) 

45900000 

001254 


I F ( IPATH.EQ .2 ) GO TO 32 

46000000 

001256 


IF( IKY.EQ.O ) PRINT 64, ( FGAM ( I ) , 1 = 1 , 15 ) , ASUM 

46100000 

001273 


64 FORMAT ( 11X 32HGAMMA CORRECTED TO MAKE ASUM=1.0/ 3 ( E28 . 8 ,4E1 7. 8/ ) 

46200000 



1 11X 5HASUM=E17 .8/ ) 

46300000 


c 


46400000 


CO 
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C 


*** ADD WATER 



46500000 

001273 



F0LD= F 



46600000 

001275 



RHAT=SQRT(VAR( 2) ) 



46700000 


C 


EQUATION < 2) 



46800000 

001277 



F=1.0 -RHAT**3 



46900000 

001301 



FWASK= F*WASK 



47000000 


C 


EQUATION (38) 



47100000 

001303 



FGAM( 1 )=FGAM( 1 ) + ( WASK/ 18 . ) * ( F-FOLD 

) 


47200000 

001307 



IF(IKY.EQ.O) PRINT 35, FGAM(l) 



47300000 

001316 


35 

FORMAT ( 1 IX 8HFGAM( 1 ) = E17.8) 



47400000 

001316 



IPATH = 2 



47500000 

001317 



GO TO 30 



47600000 

001320 


32 

CONTINUE 



47700000 

001320 



FGAM( 15)=Q. 



47800000 

001321 



DO 34 1=10,14 



47900000 

001323 


34 

FGAM ( 1 5 ) =F G AM ( 15)+FGAM( I ) 



48000000 


C 


*** PRINT GAMMAS AFTER ADJUSTMENT 

AND WATER ADDITION 

48100000 

001327 



I F ( IKY. EQ. 0 ) PRINT 33, 

( FGAM ( I ) , 1=1 , 15 ) , ASUM 

48200000 

001343 


33 

FORMAT ( 1 IX 53HGAMMA, COMPOSITION AFTER 

WATER ADDITION AND CORRECTIO 

48300000 




IN/ 3(E28.8,4E17.8/ ) 11X 5H ASUM=E17. 8/ ) 

48400000 


C 


PUT CORRECTED GAMMAS 

IN 

VAR 

48500000 


C 


OR REPLACE INTEGRATED 

GAMMAS WITH ADJUSTED, WATERED 

48600000 


c 


DOWN GAMMAS 



48700000 

001343 



DO 38 1=1,15 



48800000 

001345 



V AR ( I +3 ) =0 . 



48900000 

001347 


38 

VAR ( 1+3 ) =FG AM( I) 



49000000 

001354 



I KY= 1 



49100000 


c 


*** AT INITIAL CONDITIONS, CONTINUE 

INTEGRATION 

49200000 

001355 



IF( VAR(l).EQ.INITX) GO TO 102 



49300000 


c 


*4* IS THE CASE FINISHED 



49400000 

001357 



IF l VAR ( 1 ) . L T. XTERM) GO TO 101 



49500000 

001361 



l F ( IFIN.NE.O) GO TO 82 



49600000 


c 


TO PRINT FINAL ANSWERS 



49700000 

001362 



IPFK= IPF 



49800000 

001364 



XPRINT=VAR (1) 



49900000 

001366 



IF I N= 1 



50000000 

001367 



GO TO 57 



50100000 


c 





50200000 

001367 


82 

CONTINUE 



50300000 

001367 



PRINT 50 



50400000 

001373 


50 

FORMAT ( 10H END CASE/) 



50500000 

001373 



GO TO 1 



50600000 

001374 



END 



50700000 


APPENDIX 


000002 

000002 

000002 

000002 

000002 

000002 

000002 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


000002 

000002 

000002 


000002 

C00002 

000002 

000002 


C 

c 

c 

c 

c 


000002 

000002 

000002 

000002 

000002 

000002 

000002 

000002 


SUBROUTINE BASIC 

SUBROUTINE BASIC EVALUATES DERIVATIVES OF R HAT SQUARED, V, AND 
GAMMA 1 THRU GAMMA 14 

IBUG= 0, NO DEBUG PRINT OUTS 

NOT=0 TO GET DEBUG PRINT OUTS 
CUVAR(l) = X 

CUVAR 1 2 ) = RHATS DERI 1 ) = DRHATD DERIV R HAT SQUARED 

CUVARI3) = V DERI 2 ) = VDOT 

CUVAR ( 4 ) = GAMU) DERI3! = GAMOOT(l) 


COMMON VAR ( 20) ,DERI20) , ELE1 ( 20 ) , ELE2 l 20 ) , C I , SP EC ,NDE, ELT , NERR 
COMMON RHOL ,HL ,DN» EL, TW , CDC ,CDFM , FRO , CP VW , MASK ,R , H ASK 

COMMON T 

COMMON E1,E2,P2»U2*RH02,V2 

COMMON F,NE, Z,XI( 15) , FAC 1 , FAC2 ,U ,RHOU , RH02U2 , RHO 
COMMON KLCI 45) ,KUC l 40 ) , MB AR , FE ,CAPU, L AMD , P 

COMMON AA, AB, AC, AD, AE, AF, AG, AH, A I , A J , AKK , AL , AM , AN , AO, AP , AQ, AR, 

1 AS,AT,AU,AV, AW, AX, AY , AZ , B A , BB , BC , BD, BE , BF, BG, BH, B I , BJ , BKK , 

2 BL,BM,BN*BO,BP*BQ» BR, BS,BT 
COMMON NTLUP.VARIX (40) ,VARCP(40> 

EQUIVALENCE I A A, ALPDUM 1 1 ) ) 

DIMENSION A LPDUM ( 46 ) 


LABI IN MAIN, BASIC, FOFX 

COMMON/LABI/ 1 BUG,CU VAR, A1 , A2 , A3 , A4 , A5 , A6 , TLI Ml , TL I M2 ,TL IM3 
LAB2 IN MAIN, BASIC, CHECK, FOFX 
C0MM0N/LAB2/ENTHAL 

LAB3 IN MAIN , BASIC AND CHECK 
COMMON /LAB 3/ TPREV .KEY , HEPS , KAY 
LAB 4 IN MAIN AND BASIC 

COMMON /LAB4/ MW 1 1 5 ) , W (1 5 ) , GAMMA , HC ITR , B, DELTX , EM 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


MULT (40) , DLT(40) ,SIGLT(40) 

GAM ( 15) , GAMDOT ( 15 ) , FORT (15) 

CPOR ( 15 ) , G( 15 ) 

TKETAB(82),Z0MTAB(82) 

EMTABI20) ,CDCT A8( 20) , CDFMTB ( 20 ) , AYATABI 2 C ) , EN TAB ( 20 ) 
AK(31,BM3), CK( 3 ) ,DK( 3 ) , EK 1 3) 

A 7( 1 5 , 4 ) 

Al( 15,4) , A2( 15,4) ,A3( 15,4) ,A4(15,4) ,A5( 15,4) ,A6(15,4) 


50800000 

50900000 

51000000 

51100000 

51200000 

51300000 

51400000 

51500000 

51600000 

51700000 

51800000 

51900000 

52000000 

52100000 

52200000 

52300000 

52400000 

52500000 

52600000 

52700000 

52800000 

52900000 

53000000 

53100000 

53200000 

53300000 

53400000 

53500000 

53600000 

53700000 

53800000 

53900000 

54000000 

54100000 

54200000 

54300000 

54400000 

54500000 

54600000 

54700000 

54800000 

54900000 
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Ol 

CO 


C00002 

C00002 

000002 

C00002 

C00002 

000002 


C 


C 

c 


C00002 

C 

000002 


c 

c 

000002 


000002 


C00002 

C 

000002 


000002 


DIMENSION CUVAR{20 ) , IBUGI 20 ) 

REAL NLT 

REAL MULT,KUC,KLC 

REAL NU,NUFM,MUF,KF , NUC , NE , L AMD , MB AR, MW 
DOUBLE PRECISION VAR 

EQUIVALENCE ( CUVAR ( 1 ) , X ) , ( CUVAR { 2 ) , RHATS ) , 

1 (CUVAR (3) »V)»(CUVAR(4),GAM(1)} 

EQUIVALENCE ( DERI 1 ) ,DRHATD) , ( DER ( 2 I , VDOT ) , ( DER ( 3 ) , GAMDOT ( 1) ) 

DATA ( A7 ( I ) ,1 = 1,60) / -.8 81789073,3.51339873,2.28 796557, 

1 -A- .06 56627 1 , 3.91170120, 2.42672054, 4. 18000000,2.53040711 , 

2 -.460878520, 3.10953628, 2.72216325, 3.53897944,4.37796250, 

3 -1.15390000,-11.7338, 

4 2.21116037, 2. 308005 C7» 3.665085 52,-1.248739 2 3, 3.78027312, 

5 4.84331426, 2.68820982, 4.65791645,-2.96581444, 4.87313762, 

6 5.88646745, 4.65578522, 4.10796742,-1.15390000,-11.7338,30*0./ 

CONSTANTS FOR USE IN REACTION RATE FORMULAS 
DATA (MULT( I) , 1 = 1,40) / 10 . 9 57 , 14. 933 , 1 4 . 0 , 1 . 867 , 1 4 .483 , 10. 578 , 

1 .970,1.778, .968,10.435, 10.182, 8.471, .947, 1.875, 7.467, 

2 7.0, 5* 5.486E-4, 14.483, 9.333,10.182,0.966, 9.545, 10.435, 

3 0.968, 7.46 7,0.933, 0.941, 10.435, 10.182,9. 545,9. 545 ,C. 96 8, 10. 435, 

4 0.968,10.182,1.875 / 

DATA 1DLT( I ), 1=1,40) / 217.10, 217.34, 414.50,190.34,275.74, 

1 186.09, 31.440, 4.249,89.646, 58.206,138.93,31.015, 26.766, 

2 62.380, 117.69, 247.27, 661.94, 393.00, 617.75, 578.24,577.81, 

3 268.94, 44.186, 83.698, 84.122, 224.75, 185.24, 184.82, 39.512, 

4 39.937, 0.425, 282.96, 94.744, 46.310, 130.01, 314.40, 99.843, 

5 129.16, 178.44, 287.63 / 

DATA! SIGLT( I) , 1 = 1,40) / 6* 10.0E-16, 4*1.0E-16, 2.69E-16, 

1 3*1 . 0E-16, 0. C18E-16, 0.066E-16, 5*0.879E-16, 19* 20.0E-16/ 

DATA (TKETABt I ) , 1=1 , 8 2 ) / . 30 , .3 5 , . 40, .45 , . 50 , . 5 5 , . 60 , .65 , . 70 , . 7 5 , 

1 .80,. 85, .90, .95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 

2 1.45, 1.50, 1.5 5, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2. 00, 

3 2. 10, 2. 20, 2. 30, 2. 40,2.50, 2. 60, 2. 70, 2. 80, 2. 90, 3. 00, 3. 10, 3. 20, 

4 3.30,3.40, 3. 50,3.60,3.70,3. 80,3.90,4.00,4. 10,4.20 ,4.30,4.40, 

5 4.50,4.60, 4. 7 0,4.80, 4. 90, 5. 00,6.00, 7. 00, 8. 00, 9. 00, 10. 0,20.0, 

6 30.0,40.0, 50. 0,60.0,70.0,80.0,90. 0,100. ,200. ,300. ,400./ 

DATA ( ZOMT AB( I ),I=1,82)/ 


55000000 

55100000 

55200000 

55300000 

55400000 

55500000 

55600000 

55700000 

55800000 

55900000 

56000000 

56100000 

56200000 

56300000 

56400000 

56500000 

56600000 

56700000 

56800000 

56900000 

57000000 

57100000 

57200000 

57300000 

57400000 

57500000 

57600000 

57700000 

57800000 

57900000 

58000000 

58100000 

58200000 

58300000 

58400000 

58500000 

58600000 

58700000 

58800000 

58900000 

59000000 

59100000 

59200000 


> 

Oj 




iiiii ii mini i hi i ii mi i 



000002 

000002 

000002 

000002 


000002 

000002 

C 

000002 

c 

c 

000002 

000004 

C00006 

000011 

C 

000012 

000016 

c 

000020 

c 

c 

000023 

000043 

C 

000045 

000062 

000063 


1 2.785, 2.628,2.492,2.368,2.257,2.156,2.065,1.982,1.908, 1.841, 

2 1.780, 1.725,1 .675,1.629,1.587,1. 549, 1. 514, 1 .482, 1 .452 , 1 .424, 

3 1.399,1.375,1.353,1.333,1.314,1.296,1.279, 1.264,1.248,1.234, 

4 1.221,1.209,1.197, 1.186,1.175,1.156,1.138,1.122,1.107,1.093, 

5 1.081,1.069,1 .058, 1.048,1.039,1.030,1.022,1.014,1.007, .9999, 

6 .9932, .9870, .9811, .9755, .9700, .9649, .9600, .9553, .9507, .9464, 

7 .9422,-9382, . 9343, .9305, .9 269, . 8963 ,. 8727, . 85 38, . 8379 ,. 8242, 

8 .7432, .7005, .6718, .65 04,. 6 335,. 61 94, . 6076, . 5973 ,. 5882 ,. 532 0, 

9 .5016,-4811/ 

DATA (EMTAB(I) ,1 = 1,20)/ . 5 , . 6 , .7 , . 8, . 9, 1 . 0 , 1 . 2 , 1 . 4 , 1 . 6 , 1 . 8 ,2 . 0, 

1 2. 4, 2. 8, 3. 2, 4. 0, 5. C, 6. 0,8. 0,10.0, 12.0/ 

DATA (CDCTAB(I) ,1 = 1 ,20)/ 0. 520, . 55 1 , . 586 , . 625 , . 666 , .7 1 2 , . 80 1, 

1 .880, .929, .955, .971, .981, .969, .949, .919, . 9 10 ,. 910 , .9 10 , .910, .910/ 
DATA (CDFMTBl I ) ,1=1 ,20)/ 7.80,6.50,5.57,4.92,4.45,4.10,3.60, 

1 3. 23, 2. 98, 2. 80, 2. 68, 2. 48, 2. 36, 2. 28, 2. 17, 2. 10, 2. 05, 2. 02, 2. 0,2.0/ 
DATA (AYATAB(I) ,1 = 1,20)/ 0 . 315 , .240, . 1 82 , . 141 , . 110 , .090 , . 065, 

1 .055, .049, .047, .046, .0455, .0455, .0453, .0452, . 045 1 ,. 0449 ,. 0448 , 

2 .0447,-0447/ 

DATA (ENTABII) ,1 = 1,201/ 0.4 10 , .460, . 500 , . 545 , . 590 , . 620, .670 , . 690, 

1 .710, .715, .720,. 725,. 725,. 730, .730, .735, .735, .740 ,.745,. 745 / 

DATA PI , A,TZ,REMD/3.1415926536,6.02486E23,273. 1 6, 8 .31 696E7/ 

EXTERNAL FOFX 

KAY SET = 1 IN BASIC WHEN T.LT.O. 

IF! KAY.NE.O ) RETURN 
IF( WASK.EQ. 0. ) RHAT S = 0 . 

RH AT = SORT I RHATS ) 

MTLUP=1 

KNOWN PRESSURE DISTRIBUTION 
CALL FTLUP( X,P0P2,MTLUP ,NTLUP , VAR I X , VARDP ) 

P= POP2 * P 2 

COMPUTE F IN MAIN PROGRAM 
TEM= 1 . +F* W ASK 

5/13/67 TREAT THE PRODUCT RHO*U AS CONSTANT OVER AN INTERVAL 
EQUATION ( 6) 

U= ( P2— P 1/(2 .*RHOU )+(l./TEM)*((P2-P)/(2 . *RH02U2 1 +U2+ 

1 (V2-(1.-F)*V)*WASK) 

IFIWASK.EQ. 0. 1 V=U 

EQUATION 1 8) 

ENTH AL= (HASK+HL*F*WASK-( 1 .-F 1 /2 . *V** 2*WA SK 1 /TEM - U**2/2. 

G AM ( 1 5 ) =0. 

DO 10 1= 10,14 


59300000 
59400000 
59500000 
59600000 
59700000 
59800000 
59900000 
60000000 
60100000 
60200000 
60300000 
60400000 
60500000 
60600000 
60700000 
60800000 
60900000 
61000000 
61100000 
61200000 > 

61300000 2 

61400000 M 

61500000 g 

61600000 O 

61700000 R 

61800000 
61900000 
62000000 
62100000 
62200000 
62300000 
62400000 
62500000 
62600000 
62700000 
62800000 
62900000 
63000000 
63100000 
63200000 
63300000 
63400000 
63500000 


cn 

co 



Cn 

4 * 


000064 


10 

GAM( 15)= GAM( 15)+ GAM( I ) 

63600000 


C 



63700000 


C 


ITERATE FOR T LF ITR1 C5 

63800000 


C 



63900000 

000071 



DELTT= . 001 

64000000 

000073 



El=. 0000001 

64100000 

000074 



E2=.0OO1 

64200000 

000076 



MAXI =100 

64300000 


C 


FOF X VERY SENSITIVE TO SMALL CHANGES IN T 

64400000 

000077 



CALL ITRHT ,D ELTT , FOFX , El , E2 ,MAX I , ICODE) 

64500000 

000105 



IF ( ICODE. EQ.O J GO TO 12 

64600000 

000106 



PRINT 11, ICODE, T.DELTT, MAXI 

64700000 

000122 


11 

FORMAT ( 22H IN BAS I C- ITR1-I CGDE= 12 ,4H T=E17.8,8H DELTT=E17. 8 , 

64800000 



1 7H MAXI= I 3/ ) 

64900000 


C 


*** STOP 11 

65000000 

000122 



STOP 11 

65100000 


C 



65200000 

000124 


12 

CONTINUE 

65300000 

000124 



IFIT.GT.O. ) GO TO 13 

65400000 

000127 



KA Y= 1 

65500000 

000130 



RETURN 

65600000 

000131 


13 

CONTINUE 

65700000 

000131 



SORTT = SQRT ( T) 

65800000 

000134 



UMV=U-V 

65900000 

000136 



TEM5= UM V*UM V 

66000000 

000137 



GAMMA=0. 

66100000 

000140 



DO 15 1 = 1,15 

66200000 

000141 


15 

GAMMA= GAMMA + GAM ( I) 

66300000 


C 



66400000 


C 


***** COMPUTE RHO, { F SUB I)X(RT), CAPITAL Kill, 1 = 1, 40, 

66500000 


c 



66600000 


c 


EQUATION (11) 

66700000 

000145 



RHO=P/ (REMD*T * GAMMA ) 

66800000 


c 



66900000 

000151 



IF { I BUG ( 1) *EQ. 0)G0 TO 17 

67000000 


c 



67100000 


c 



67200000 

000152 



PRINT 18, P, POP2, T 

67300000 

000163 


18 

FORMAT < * P , POP2 »T = * 3E17.8) 

67400000 

000163 



PRINT 16,F, TEM,U,ENTHAL, GAMMA, RHO, IGAM (I) , 1 = 1 , 1 5 ) 

67500000 

000211 


16 

FORMAT ( 26H F , TEM , U , ENTH AL ,GAMM A , RH0/6E 17 . 8 /6H GAM=/3 ( 5E20 .8/) ) 

67600000 


c 



67700000 

000211 


17 

CONTINUE 

67800000 
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000211 


NLT = AL3G(T) 

67900000 

000214 


T2= T**2 

68000000 

000215 


T3= T2*T 

68100000 

000216 


T4= T3*T 

68200000 


C 

TLIM1=1000. 

68300000 


c 

TL IM2=1 500 0 . 

68400000 


c 

TL IM3=1000000. 

68500000 

000217 


IF( T . LE .TL I Ml ) IT = 1 

68600000 

000222 


IF(T.GT.TLI Ml .AND.T.LE.TLIM2) IT = 2 

68700000 

000234 


IF(T .GT .TLIM2.ANO.T .LE.TLIM3) IT=3 

68800000 

000246 


IF(T.GT.TLIM3) IT=4 

68900000 

000251 

c 

DO 20 1=1, 15 

69000000 

69100000 


c 

(F SUB I)/RT COEFFICIENTS FROM ZELEZNIK AND GORDON 

69200000 


c 

c 

NASA TN D- 1454 

69300000 

69400000 

000253 


20 FORT ( I ) = A1(I,IT)*(1.-NLT)-A2(I,IT)/2.*T- A3 ( I , IT) /6. *T2-A4 ( I , I T) 
1 / 1 2. *T3 - A 5( I * IT ) / 20 . *T4 +A6(I,IT)/T -A7(I,IT) 

69500000 

69600000 


c 

FIND CAP K (KUC) 

69700000 


c 

R, UNIVERSAL GAS CONSTANT, 82.1023(ATM CM**3)/ DEG. K 

69800000 


c 

R, UNIVERSAL GAS CONSTANT, 8.31696E7 ERG/ ( MOLE-DEG.K) 

69900000 

000322 


OORT= 1 . / t 82. 1023*T) 

70000000 

000325 


KUC(1)= EXP (— FORT ( 6 )-FORT( 9 ) +FORT ( 1 ) ) *OORT 

70100000 


c 

EQUATION (231 

70200000 

000334 


KUC ( 2 )= EXP (— F ORT (8)-F0RT(8) + F0RT(2) ) *OORT 

70300000 

000342 


KUC ( 3 )= EXP (— F0RT(7)— FORT(7J+FORT( 3) ) *OORT 

70400000 

000350 


KUC C 4 ) = EXP (-FORT ( 9 I-FORT ( 9 ) +FORT (4) ) *00RT 

70500000 

000356 


KUC ( 5 )= EXP (-FORT(7)-FORT(8)+FORT( 5) ) *OORT 

70600000 

000365 

c 

KUC ( 6 ) = EXP (-FORT ( 8 I-FORT (9 ) +FORT ( 6) ) *OORT 

70700000 

70800000 

000374 


KUC (17)=EXP (-FORT ( 11)-F0RT( 15)+ F0RT(3 II *OORT 

70900000 

000403 


KUC(18)=EXP(— FORT! 10)-F0RT( 151+ F0RT(5 )) *OORT 

71000000 

000412 


KUC ( 1 9 >=EXP (-FORT ( 1 2 1 -FORT ( 1 5 ) + F0RTI7 )) *OORT 

71100000 

000421 


KUC(20)=EXP (-FORT! 13 ) -FORT (15)+ F0RT(8 )) *OORT 

71200000 

000430 

c 

KUC(21)=EXP (— FORT( 14)-F0RT( 15)+ F0RT(9 )) *OORT 

71300000 

71400000 

000437 


KUC ( 7 ) = E XP ( - FORT ( 6 ) -FORT(8 ) +FORT(2 ) +F0RT(9 )) 

71500000 

000446 


KUC ( 8 ) = EXP(— FORT ( 6 ) -FCRT(9 ) +F0RT(4 ) +FORT(8 )) 

71600000 

000455 


KUC ( 9 ) = EX P ( - FORT{ 6 ) -FORT(7 ) +FORTI 5 ) +F0RT19 )) 

71700000 

000464 


KUC (10)= E X P( - FORT ( 2 ) -F0RT(7 ) +FORT(5 ) +F0RT(8 )) 

71800000 

000473 


KUC (11)= EX P ( - FORT ( 5 ) -F0RT(7 ) +FORT(3 ) +FORT( 8 )) 

71900000 

000502 


KUC (12)= EX P ( - FORT ( 6 ) -F0RT(6 ) +FORT( 1 ) +FORT(8 )) 

72000000 

000510 


KUC( 13 )= E X P (— FORT ( 4 ) -F0RT(6 ) +FORT11 ) +F0RTI9 )) 

72100000 


cn 

cn 
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Ol 

Ol 


000517 

000526 

000535 

000543 

000552 

000561 

000570 

000577 

000606 

000615 

000624 

000633 

000642 

000651 

000660 

000667 

000676 

000705 

000714 

000723 

000732 

000741 

000750 

000751 

000776 


000776 
000776 
001000 
CO 100 3 
001005 
001012 
001015 
001020 

001021 

001026 


C 


C 

C 

C 

C 


KUC (14)= EXP ( -FORT ( 1 ) -F0RT(7 ) +F0RT(4 ) +F0RT(5 )) 

KUC (15)= E X P( -FORT 110) -F0RTI15) +F0RT<7 ) +F0RK8 )) 

KUC (16)= EXP(-FORHll) -F0RK15) +F0RT(7 ) +F0RK7 )) 


KUC ( 22 ) = E X P ( - FORT ( 5 ) -FORT(ll) +F0RT(3 ) + FORT(IO)) 

KUC (23)= E X P( - FORT { 7 ) -FORT(ll) +F0RT13 ) + F0RT112)) 

KUC ( 2 4 ) = E X P ( - FORT { 8 ) -FORT (11) +F0RTI3 ) + F0RK13)) 

K UC ( 2 5 ) = E X P ( - FORT ( 9 ) -FORT (11) +F0RT(3 1 + F0RT(14)) 

KUC ( 26 ) = E X P ( - FORT ( 5 ) -F0RT112) +F0RTI7 ) + FORT(IO)) 

KUC ( 27 ) = E X P ( - FORT ( 5 ) -F0RT(13) +F0RT{8 ) + FORT(IO)) 

KUC ( 2 8 ) = E X P ( - FORT ( 5 ) -F0RTI14) +F0RT(9 ) + FORT(IO)) 

KUC ( 2 9)= E X P ( - FORT ( 8 ) -F0RT(12) +F0RTI7 ) + F0RT(13)) 

KUC ( 3 0 ) = E X P ( - FORT ( 9 ) -F0RT(12) +F0RT(7 ) + F0RT(14)) 

KUC (31)= EXP ( -FORT ( 9 ) -F0RTI13) +F0RT<8 ) + F0RT(14)> 
KUC ( 32 ) = E X P ( — FORT ( 2 ) -F0RT(12) +F0RT(8 ) + FORT(IO)) 

KUC ( 3 3 ) = E X P (- FORT ( 5 ) -F0RT(12) +FORT ( 8 ) + FORT (11)) 

KUC ( 34) = E X P ( - FORT( 3 ) -F0RT(13) +F0RT(7 ) + FORT(IO)) 
KUC ( 3 5 ) = E X P ( - FORT ( 8 ) -FORT(ll) +F0RT(7 ) + FORT(IO)) 

KUC ( 36 ) = E X P ( - FORT ( 6 ) -F0RT(12) +F0RTI9 ) + FORT(IC)) 

KUC ( 3 7 ) = EX P ( - FORT ( 2 ! -F0RT(12) +F0RT(5 ) + F0RT(13)) 

KUC (38)= EXP(-F0RT{6 ) -F0RT(12) +F0RT15 ) + F0RT(14)( 

KUC ( 39 ) = E X P ( - FORT ( 5 ) -F0RT(12) +F0RT<3 ) + F0RT(13)) 


KUC ( 40 ) = E X P ( - FORT ( 1 ) -F0RTI12) +F0RT(4 ) + FORT(IO) 


IF( IBUG(2) • EO.O) GO TO 22 

PRINT 21,OORT,NLT,T, ( KUC ( I ) , I =1 , 40 ) , ( FORT ( I ) , I =1 , 1 5 ) 

21 FORMAT ( 13T OORT ,NL T , T=3E 20 . 8/ 6H KUC=/8 < 5E20 . 8/ ) , 7H FORT= / 
1 3 ( 5E20 .8/ ) ) 


FIND KLC, SPECIFIC REACTION RATE CONSTANTS 
A=AVOGA DRO CONSTANT PARTICLES/MOLE 
LOWER CASE K( L) ,L=1 ,40 

22 CONTINUE 
TZOT =T Z/T 

EREMDT = 8.0*REM0*T/P I 

DO 30 1= 1,40 

F AC3= SORT ( EREMOT ZMULT ( I ) ) 

DLTTZOT- DLT ( I ) *TZOT 
FAC4=EXP(-DLTTZ0T) 

FAC34=FAC3*FAC4 

EQUATION (18) 

32 KLC( I )=A*SI GLT ( I )*FAC34 
IF( I .GT.l) GO TO 31 


72200000 

72300000 

72400000 

72500000 

72600000 

72700000 

72800000 

72900000 

73000000 

73100000 

73200000 

73300000 

73400000 

73500000 

73600000 

73700000 

73800000 

73900000 

74000000 

74100000 

74200000 

74300000 

74400000 

74500000 

74600000 

74700000 

74800000 

74900000 

75000000 

75100000 

75200000 

75300000 

75400000 

75500000 

75600000 

75700000 

75800000 

75900000 

76000000 

76100000 

76200000 

76300000 

76400000 
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C 


EQUATION (16) 

76500000 

001030 



KLCII ) = KLC m*(FAC2*(DLTTZOT + 1.0)**2) 

76600000 

001035 



GO TO 34 

76700000 

001035 


31 

IF( I .GT.6) GO TO 34 

76800000 


c 


EQUATION (17) 

76900000 

001041 



KLC ( I )= KLC (I )*(FAC1*(DLTTZ0T+0.5J**1.5) 

77000000 

001050 


34 

I F ( I BUG ( 31 . EQ.O) GO TO 30 

77100000 

001051 



PRINT 33, FAC1,FAC2,FAC3,FAC4,SIGLT(I),DLT( I),MULT(I) 

77200000 

001076 


33 

FORMAT ( 7E17 .8 ) 

77300000 

001076 


30 

CONTINUE 

77400000 


c 



77500000 


c 


***** COMPUTE M BAR, N SUB E, LAMBDA SUB D, CAP U, 

77600000 


c 


LOWER CASE K ( L ) ,L=41,45 

77700000 


c 



77800000 

001100 



MBAR=0 

77900000 

001101 



DO 40 I =10,14 

78000000 

001103 


40 

MBAR=MBAR + GAM(I)* MW ( I ) 

78100000 


c 


EQUATION (31) 

78200000 

001111 



MBAR= MBAR/ GAM(15) 

78300000 


c 


EQUATION (40) 

78400000 

001112 



NE=A*RHO*GAM( 15) 

78500000 


c 


EQUATION (30) 

78600000 

001115 



LAMD= 6.90* SQR T( T/N E) 

78700000 


c 


EQUATION (29) 

78800000 

001122 



CAPU= UMV/( SQRT(5.486E-4/MBAR)*6.22E5*SQRTT ) 

78900000 

001131 



IFtWASK.EQ. 0. ) VTEST=0 . 

79000000 

001133 



IF(WASK.EQ. 0. ) GO TO 35 

79100000 

001134 



VTE ST= ( 600 . *RH02U2*WASKI/( PI*RHAT**3*RH0L) 

79200000 

001142 


35 

CONTINUE 

79300000 

001142 



I F ( V.LE.VTEST) GO TO 43 

79400000 

001145 



FIRST=(3.*U#RHATS) / ( 4 . *R*RHOL*V*TE M ) *WASK 

79500000 


c 


EQUATION (255 

79600000 

001153 



FIRST=FIRST* 6.22E5*SQRTT*FE 

79700000 

001157 



DO 45 1=41 , 45 

79800000 

001160 


45 

KLC ( I ) = FIR S T 

79900000 

001164 



GO TO 42 

80000000 

001165 


43 

DO 44 1=41,45 

80100000 

001167 


44 

KLC ( I )=0. 

80200000 


c 



80300000 


c 



80400000 

001172 


42 

CONTINUE 

80500000 

001172 



IF ( I BUG ( 3) . EQ. 0)G0 TO 47 

80600000 

001173 



PRINT 46, (KLC(I) ,1=1,45), MBAR, NE,L AM C,CAPU,GA,GB,GL,FE,VTEST 

80700000 


cn 

-a 



U1 

oo 


001227 46 FORMAT ( 6H KLC=/9 ( 5 E20 . 8/ ) / 20H MB AR , NE , L AMD , C A PU= 4E2 0 . 8/19H GA,G 

1B,G,FE,VTEST=5E18.8/) 

001227 PRINT 41 »R » RH AT 

001237 41 FORMATI2E17 .8) 

C 

C AA THRU BS 

001237 47 CONTINUE 

001237 AA= KLCll)*(RH0*GAM(6)*GAM(9)*GAMMA/KUC( 1) -GAM ( 1 ) *GAMMAJ 

C 

C AB=KLC( 2)*<B2,2-G2,2) 

C B2,2= EQUATION (19) 

C G2, 2= EQUATION (20) 

001246 AB= KLC ( 2 ) * ( ( RHO *G AM ( 8 ) *GA M ( 8 ) * GAMM A ) /KUC ( 2 ) - GAM( 2 ) *GAMMA ) 

001254 AC= KLC ( 3 ) * ( ( RHO *G AM ( 7 )*GAM ( 7 ) * G AMMA ) /KUC ( 3 ) - GAM( 3 ) *GAMMA ) 

001262 AD= KLC ( 4) * ( ( RHO *G AM (9 ) *GAM ( 9 ) * GAMM A ) /KUC ( 4 1 - G AM ( 4 ) *GAMMA ) 

001270 AE= KLC(5)*((RH0 *G AM( 7 ) *GAM ( 8 ) * GAMMA ) /KUC ( 5 ) - GAM( 5 ) *GAMMA ) 

001276 AF= KLC (6) * ( ( RHO *G AM ( 8 ) *GAM ( 9 ) * GAMMA ) /KUC ( 6) - G AM ( 6 ) *GAMMA ) 

001304 AG= KLC ( 7) * U GAM ( 6 )* G AM ( 8 ))/KUC(7 ) - GAM ( 2 ) *GAM(9 ) ) 

001311 AH= KLC ( 8 )*( ( GAM ( 6 )* G AM ( 9 ))/KUC(8 ) - GAM ( 4 ) *GAM(8 ) ) 

001316 A 1= KLC ( 9 )*{( GAM ( 6 )* GAM ( 7 ))/KUC(9 ) - GAM ( 5 ) *GAM(9 ) ) 

001324 AJ= KLC ( 10 ) *( ( GAM ( 2 )* GAM (7 ))/KUC(10) - GAM( 5 ) *GAM(8 I ) 

001331 AKK=KLC( 11 ) *( ( GAM ( 5 )* GAM (7 ))/KUC(ll) - GAM ( 3 ) *GAM<8 ) ) 

001336 AL= KLC ( 1 2 )* ( (GAM ( 6 )*GAM(6 ))/KUC(12) - GAM( 1 )* GAM( 8 ) ) 

001343 AM= KLC(13 )*( (GAM( 4 )*GAM(6 ))/KUC(13) - GAM( 1 )* GAM ( 9 ) ) 

001350 AN= KLC ( 1 4 )* ( ( GAM ( 1 )*GAM(7 ) ) /KUC 1 14) - GAM( 4 )* GAM ( 5 ) ) 

001355 A0= KLC ( 15 )* ( (GAM (10) *GAM (151) /KUC (15) - GAM ( 7 )* GAM ( 8 ) ) 

C 

C AP=KLC( 16 ) *{ B7» 16-G7»16) 

C 67,16= EQUATION (21) 

C G7 ,16= EQUATION (22) 

001363 AP= KLC ( 1 6 ) * ( (GAM( 11 ) *GAM ( 15) ) /KUC(16) - GAM( 7 )* GAM( 7 ) ) 

001370 AQ=KLC( 17) * ( (RHO*GAM( 11)*GAM( 15)*GAM( 15) ) /KUC ( 17 ) -GAM ( 3 )*GAM( 15) ) 

001376 AR=KLC( 18) * ( ( R HO*G AM ( 10 1 *GAM ( 1 5 ) *GAM ( 15 ) ) /KUC( 18) -GAM ( 5 )*GAM( 15) ) 

001404 A$ = KLC(19)*((R HO*GAM( 12)*GAM(15)*GAM(15) ) /KUC ( 19) -GAM ( 7 )*GAM( 15)) 

001412 AT=KLC(20)*((RH0*GAM(13)*GAM( 15)*GAM( 15 ) ) /KUC ( 20) -GAM ( 8 )*GAM( 15) ) 

001420 AU=KLC(21)* ( (RHO*GAM( 14 ) *GAM ( 15 ) *GAM( 15) ) /KUC ( 21 ) -GAM ( 9 )*GAM( 15) ) 

001426 AV=KLC(22)*((GAM(5 ) *G AM ( 1 1 ) ) /KUC ( 22 ) - GAM( 3 )*GAM(10) ) 

001433 AW=KLC(23)*((GAM(7 ) *GAM ( 11 ) ) /KUC ( 23 ) - GAM( 3 )*GAM(12) ) 

001441 AX=KLC(24)*((GAM(8 ) *GAM( 11 ) ) /KUC ( 24) - GAM(3 )*GAM(13) ) 

001446 AY=KLC( 25) * ( (GAM(9 )*GAM( 11 )) /KUC (25) - GAM( 3 )*GAM(14) ) 

001454 AZ=KLC ( 26) * ( ( GAM ( 5 ) *GAM ( 1 2 ) ) /KUC ( 26 ) - GAM(7 )*GAM(10) ) 

001461 BA=KLC(27)4((GAM(5 )*GAM ( 13 ) ) /KUC ( 27 ) - GAM( 8 )*GAM(10) ) 

001467 BB=KLC (28)* ( (GAM (5 ) *GAM ( 14 ) ) /KUC ( 2 8 ) - GAM(9 )*GAM(10) ) 


80800000 

80900000 

81000000 

81100000 

81200000 

81300000 

81400000 

81500000 

81600000 

81900001 

81900002 

81900003 

82000000 

82100000 

82200000 

82300000 

82400000 

82500000 

82600000 

82700000 

82800000 

82900000 

83000000 

83100000 

83200000 

83300000 

83400000 

83700001 

83700002 

83700003 

83800000 

83900000 

84000000 

84100000 

84200000 

84300000 

84400000 

84500000 

84600000 

84700000 

84800000 

84900000 

85000000 


> 


K 



001474 


BC=KLC(29)*((GAM(8 ) *G AM ( 12 ) ) /KUC C 29 ) - GAM( 7 )*GAM(13) ) 

85100000 

001501 


BD=KLC(30)*((GAM<9 ) *G AM ( 12 ) ) /KUC ( 30 ) - GAM( 7 )*GAM(14) ) 

85200000 

001507 


BE=KLC(31)*((GAM(9 ) *GAM ( 13 J ) /KUC ( 3 U - G AM( 8 J*GAM(14) ) 

85300000 

001514 


BF=KLC(32)*((GAM(2 ) ♦GAM ( 12 ) 1 /KUC ( 32 ) - GAM( 8 )*GAM(10) ) 

85400000 

001521 


BG=KLC ( 33) ♦ { ( G AM ( 5 )*GAM( 12 ) ) /KUC ( 33 ) -GAM18 )*GAM(11) ) 

85500000 

001527 


BH= KLC ( 34 ) *( ( GAM ( 3 ) *GAM ( 1 3 ) ) /KUC { 34 )-GAM ( 7 )*GAM(10) ) 

85600000 

001534 


BI= KLC ( 35 ) *( ( GAM( 8 ) *GAM ( 1 1 ) ) /KUC ( 35 )-GAM( 7 )*GAM{10) ) 

85700000 

001541 


BJ= KLC(36) *{ <GAM(6 ) *GAM( 12 ) J/KUC ( 36 )-GAM( 9 )*GAM{10) ) 

85800000 

001546 


BKK=KLC ( 37) *t ( GAM ( 2 ) *G AM { 12 ) ) /KUC ( 37 )-GAM ( 5 )*GAM(13) ) 

85900000 

001553 


BL= KLC ( 38 ) *{ ( GAM( 6 ) *GAM( 12 ) J/KUC { 38 )-GAM( 5 )*GAM(14) ) 

86000000 

001561 


BM= KLC ( 39 ) *{ ( GAM( 5 ) *GAM ( 1 2 ) ) /KUC { 39 ) -GAM ( 3 )*GAM(13) ) 

86100000 

001566 


BN= KLC (40 ) ♦( ( GAM( 1 )*GAM(12J ) /KUC ( 40 )-GAM( 4 >*GAM(10> ) 

86200000 

001574 


B0= KLC (41) ♦ GAM(IO) 

86300000 

001575 


BP= KLC (42) ♦ GAM(ll) 

86400000 

001577 


B0= KLC (43) ♦ GAM (12) 

86500000 

001601 


BR= KLC (44) * GAM( 13) 

86600000 

001603 


BS= KLC (45 ) * GAM( 14) 

86700000 


C 


86800000 

001605 


I F ( IBUG(4) . EQ.OIGO TO 49 

86900000 

001606 


PRINT 48, AA,AB,AC,AD,AE,AF,AG»AH,AI, A J , AKK , AL , AM , AN , AO, AP, AQ,AR, 

87000000 



1 AS, AT,AU,AV, AW, AX, AY , AZ, BA , BB , BC , BD, BE , BF , BG, BH, B I , B J , BKK, BL,BM, 

87100000 



2 BN,BO,BP,BQ,BR,BS 

87200000 

001744 


48 FORMAT ( 20H A A , AB, AC, AD, AE , AF=/6E20. 8/20H AG , AH, A I ,A J , AK, A L=/ 

87300000 



1 6E20.8/20H AM, AN , AO , AP , AQ , AR=/6E20. 8/20H AS , AT , AU , AV , AW, AX=/ 

87400000 



2 6E20.8/20H A Y , AZ , B A, BB , BC , BD=/6E20. 8/20H BE , BF , BG, BH ,BI , BJ=/ 

87500000 



3 6E20.8/20H B K , BL , BM , BN , BO , BP=/6E20 . 8/ 1 1 H BQ , BR , BS=/ 3E20. 8/ ) 

87600000 


C 


87700000 


C 

BEGIN RHAT SQUARED DOT EVALUATION 

87800000 


C 

COMPUTE (CP/R) SUB I 

87900000 


c 

TL I Ml = 1 000 • DEG.K 

88000000 


c 

TL I M 2=15000 • DEG K 

88100000 

001744 


49 CONTINUE 

88200000 

001744 


IF( T. LE.TLI Ml ) IT=1 

88300000 

001750 


IF ( T .GT .TL I Ml . AND.T .LE.TLIM2) IT=2 

88400000 

001762 


IF(T.GT.TLIM2. AND.T .LE.TLIM3) IT=3 

88500000 

001774 


IF(T.GT.TLIM3) IT=4 

88600000 

001777 


DO 51 1 = 1,15 

88700000 


c 


88800000 


c 

CP/R 

88900000 


c 


89000000 

002001 


51 CPOR ( I ) = A1 ( I, IT) +A 2 ( I , IT ) *T + A3(I,IT)*T2 +A4(I,IT)*T3 

89100000 



1 + A5(I , IT) *T4 

89200000 

002031 


CPBAR= 0. 

89300000 


U1 

CD 


APPENDIX 



05 

O 


002032 



00 50 1=1,15 

89400000 

002033 


50 

CPB AR= CPBAR+ GAM(I)* CPdRU) 

89500000 


C 


EQUATION (50) 

89600000 

002041 



CPB AR= CPB AR *REMD 

89700000 

002042 



TS= T+TEM5/ (2.*CPBAR) 

89800000 


c 


EQUATION (49) 

89900000 

002046 



T AV = ( TS+TW) /2. 

90000000 

002050 



TAV809= TAV/809.1 

90100000 

002052 



CALL FTLUP ( TAV809, OMEGA, 1,82,TKETAB, ZOMTAB) 

90200000 


c 


FIND OMEGA 2.641**2=6.974881 

90300000 


c 


EQUATION (48) 

90400000 

002056 



MUF = ( 26. 693*SQRT{ 18.*TAV) )/ ( 6.974881*0MEGA ) *1 . OE-6 

90500000 


c 


EQUATION (45) 

90600000 

002067 



RES=(RHO*UMV*l 2.*R*RHAT) ) /MUF 

90700000 


c 


EQUATION (59) 

90800000 

002074 



NUC= 2.+ .6 * SORT ( RES ) 

90900000 


c 


EQUATION (61) 

91000000 

002101 



DEL T= ( FRO*T EM5 ) / ( 2 . *CPB AR ) +T-TW 

91100000 


c 


FIND (CP/R) FOR H20 AT T=TA V R/ 1 8 .=. 462053 33E7 

91200000 

002106 



IF(TAV.LE.TLIMl) IT=1 

91300000 

002112 



IFfTAV.GT.TLIMl .AND .T AV . LE . TL I M2 ) IT=2 

91400000 

002124 



IF(TAV.GT.TLIM2 . AND .T AV . LE . TL IM3 ) IT=3 

91500000 

002136 



IF( TAV.GT.TLIM3) IT=4 

91600000 

002141 



CPOR SP= A1 ( 1 , IT )+A2( 1, IT)*TAV+A3 (1, IT ) *T AV**2+ A4 ( 1 , IT) *TAV**3 

91700000 




L +A5( 1 , I T) *TAV**4 

91800000 


r 


EQUATION (57) 

91900000 

002163 



CPF=.46205333E7*CPORSP 

92000000 

002164 



TEM= 1.+ ( C PF* DELT ) /EL 

92100000 

002170 



TEM1= ALOG(TEM) 

92200000 


c 


EQUATION (60 

92300000 

002172 



Q0GZ=EL/(CPF*DELT)*TEM1 

92400000 


c 


2*R=16.63392E7 ERG/MOLE-DEG 

92500000 


c 


EQUATION (65) 

92600000 

002175 



S=UM V/SQRT ( 16 .6 3392 E7*T*GAMMA) 

92700000 

002203 



SUMGG=0 . 

92800000 

002204 



DO 55 1= 1,15 

92900000 


c 


EQUATION (47) 

93000000 

002205 



G ( I ) = CPOR (I ) /(CPOR( I )-l .) 

93100000 

002211 


55 

SUMGG= SUMGG + GAM(I)*G(I> 

93200000 


c 


EQUATION (56) 

93300000 

002217 



KF = CPF * MUF 

93400000 

002221 



IF( S.NE.O. ) GO TO 52 

93500000 

002222 



FR=0 . 

93600000 
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002223 


ST=0. 

93700000 

002224 


NUFM=0. 

93800000 

002225 


NU=0 . 

93900000 

002226 


GO TO 53 

94000000 

002226 


52 CONTINUE 

94100000 

002226 


TEM2= EXP(-.707*S) 

94200000 


C 

EQUATION (63) 

94300000 

002233 


FR = SUMGG/ ( SUMGG+GAMMA) * (2.+.7*TEM2) 

94400000 

002240 


TEM3= EXP( -1.1 66* S**.406) 

94500000 


C 

EQUATION (64) 

94600000 

002247 


ST = . 1 041 / S** 1.14 +. 12 5* ( 1 TEM3 ) 

94700000 

002257 


ST1=S** .406 

94800000 

002263 

C 

512=5**1. 14 

94900000 

95000000 

002267 


TEM4= CPBAR *( T-TW) 

95100000 

002272 


TEM6 = ( (FR*TEM5)/2. + TEM4) /( { FRO*TEM5 ) /2 . + TEM4 ) 

95200000 

002301 


TEM7= (RHO*UMV*2.*R*RHAT*CPVW*ST)/KF 

95300000 


c 

EQUATION (62) 

95400000 

002310 


NUFM = ( SUMGG+GAMMA ) /SUMGG * TEM7 * TEM6 

95500000 


c 

EQUATION (58) 

95600000 

002314 


NU= (NUC *Q0QZ)/( 1.+ NUC*QOQZ/NUFM ) 

95700000 

002322 


53 CONTINUE 

95800000 

002322 


IF ( WASK . EQ . 0. ) DRHATO=C . 

95900000 

002324 


IF(WASK.EQ.O.) GO TO 58 

96000000 

002325 


I F ( V.LE.VTEST) DRHATD=0 . 

96100000 

002331 


IF ( V.LE.VTEST) GO TO 58 

96200000 

002333 


TEM8=(-NU*KF*0N) /( EL*R**2*RHOL*CP8AR*V) 

96300000 


c 

EQUATION (55) 

96300001 

002343 


TEM9= ( FRO*T EM 5 ) / 2 . +CP B AR*T-CPBAR*TW 

96400000 


c 

(R HAT ) **2 DOT = DER(l) 

96500000 


c 

EQUATION (54) 

96600000 

002351 


DRHATO = TEM8* TEM9 

96700000 

002352 

c 

58 CONTINUE 

96800000 

96900000 

002352 

c 

RH02=RH0**2 

97000000 

97100000 

002354 


I F ( I BUG (5) . EQ. 0 ) GO TO 57 

97200000 

002355 


PRINT 56, (CP0R( I), 1=1, 15), CPBAR, TS,TAV,MUF, OMEGA, NUC, DELT, CPF, 

1 TEM, TEM1, Q OQZ,S, SUMGG, TEM2 , FR, KF , TEM3 , ST, TEM4 , TEM5 ,T EM6 ,TE M7 , 

2 NUFM,NU»TEM8 , TEM9 , RHATDD , RHATD , DRHATO , BT ,R HO 2, RES 

97300000 

97400000 

97500000 

002466 


56 FORMAT ( 7H CPO R=/3 ( 5E20 . 8 / ) /29H C P8 AR , T S , T A V , MUF , OMEG A ,NUC=/ 

1 6E20.8/27H D EL T , C PF , TEM ,TEM1 , QOQ Z , S= / 6E 20 . 8 / 27H SUMGG, TEM2 , FR, K 
2F,TEM3,ST = /6E2 0.8/3 OH TE M4 ,TEM5 ,T EM6 ,T EM7, NUFM ,NU=/6E20.8/35H TE 

97600000 

97700000 

97800000 


a> 
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to 




3M8,TEM9,RHATDO,RHATD,DRHATD,BT=/6E20.8/7H RH02=E2C.8, 

97900000 



4 6K R E S=E 18.8/) 

98000000 


C 


98100000 


C 


98200000 


C 


98300000 


C 


98400000 


C 

COMPUTE W SUB I , 1 = 1 ,14 

98500000 

002466 


57 CONTINUE 

98600000 


c 

EQUATION (14) 

98700000 

002466 


W<1)= RH02* (AA+AL+AM-AN-BN ) 

98800000 

002474 


W ( 2 ) = RH02* (AB+AG-AJ-BF-BKK) 

98900000 

002502 


W ( 3 ) = RH02* ( AC+AK.K + AQ + AV+ AW + AX+AY-8H+BM + 8P ) 

99000000 

002515 


W ( 4 ) = RH02* (AO+AH-AM+AN+BN) 

99100000 

002523 


W I 5 ) = RH02* { AE+AI+A J-AKK+AN+AR-AV-AZ-BA-BS-BG+BKK+BL-BM+BO) 

99200000 

002543 


W ( 6 ) = RH02* ( -A A+ AF-AG- AH- AI-2 . * AL - AM-B J-8 L ) 

99300000 

002557 


W ( 7 ) = RH02*(-2.*AC-AE-AI - A J- AKK- AN+ AO+2 .*AP+AS-AW+AZ+BC+BD+BH+BI 

99400000 



1 +BQ) 

99500000 

002602 


W( 8) = RH02* ( - 2 .*AB- AE- AF- AG+AH+A J+AKK+AL+ AO+AT-AX+8 A-BC+BE+BF+BG 

99600000 



1 -BI+BR) 

99700000 

002626 


W ( 9 ) = RH02*(-AA-2.*AD-AF + AG-AH+AI+AM + AU-AY+-BB-BD-BE+BJ+BS) 

99800000 

002646 


W( 10)= R HO 2 * ( -A0-AR + AV+AZ + 8A*-BB+BF +BH + B I + B J + BN-BO 5 

99900000 

002663 


W( 11 )= RH02 * ( - AP-AQ-AV-AW- AX- AY+BG-B I-BP ) 

100000000 

C02675 


Wl 12)= RH02*( - AS+AH-AZ-BC-BC-BF-BG-BJ-BKK-BL-BM-BN-BQ) 

100100000 

002713 


W( 13) = RH0 2 *( — AT+ AX-BA+BC-BE-BH + BKK+BM-BR ) 

100200000 

002725 


W( 14)= RH02 *1 -AU+AY-BB + BD+BE+BL-BS ) 

100300000 


c 


100400000 


c 

GAMDOT ( 1 ) =DER<4) - DER(18) 

100500000 

002735 


DO 60 1 = 1 , 14 

100600000 


c 

EQUATION (13) 

100700000 

002736 


60 GAMDOT ( I ) = (DN* W ( I ) ) /(RHO*U) 

100800000 


c 

EQUATION (46) 

100900000 

002747 


EM= UMV/SQRT(REMD*T*SUMGG) 

101000000 

002755 


IF(WASK.NE.O. ) GO TO 67 

101100000 

002756 


CDZ=0. 

101200000 

002757 


C D= 0 .. 

101300000 

002760 


GO TO 65 

101400000 

002760 


67 CONTINUE 

101500000 

002760 


TEM10= EXP( -.028*RES**.82) 

101600000 


c 

EQUATION (44) 

101700000 

002767 


CDZ = 0.4 + 24. /RES +1.6+TEM10 

101800000 


c 


101900000 

002774 


68 CONTINUE 

102000000 

002774 


IF ( EM.GE..5) GO TO 61 

102100000 



002777 



FMDRES=CDZ+51. 1*1 EM /RES) 

102200000 

003002 



TEM11=1 .0+0 .256*EM*FM0RES 

102300000 


C 


EQUATION (43) 

102400000 

003006 



CD= FMORES/ TEM11 

102500000 

003007 



GO TO 65 

102600000 

003007 


61 

CALL FTLUPI EM* AY A , 1 , 20 , EMT AB , A YAT AB ) 

102700000 

003013 



CALL FTLUPI EM, EN , 1 , 2 0 , EMT AB , ENT AB ) 

102800000 

003017 



CALL FTLUPI EH, CDC , 1 , 20 , EMT AB , CDCT AB ) 

102900000 

003023 



CALL FTLUPI EM , COFM , 1 , 20 , EMT AB ,C DFMT B) 

103000000 

003027 



TEM11=EXP(-AYA*RES**EN) 

103100000 


C 


EQUATION (51) 

103200000 

003037 



CO = CDC + (CDFM -CDC)* TEM11 

-103300000 

003043 


65 

CONTINUE 

103400000 

003043 



IF! IBUG ( 6) % EQ . 0) GO TO 63 

103500000 

003044 



PRINT 62, EM, AY A, EN, CDC, CDFM 

103600000 

003062 


62 

FORMAT 1 2 1H EM , AYA , EN ,C DC ,CDFM=5E 1 8 . 8/ ) 

103700000 

003062 


63 

CONTINUE 

103800000 


C 



103900000 


C 


*** VDOT = DERI2) *** 

104000000 

003062 



IFI WASK.EQ. 0. ) VDOT = 0 . 

104100000 

003064 



IF! WASK.EQ. 0. ) GO TO 64 

104200000 


C 


DERIl) = DRHATD . DERI V R HAT SQ. 

104300000 

003065 



IFIV.LE.VTEST) V00T=0. 

104400000 

003071 



IFI V.LE.VTEST) GO TO 64 

104500000 


C 


EQUATION (42) 

104600000 

003073 



VDOT= ( 3.*DN*RH0*TEM5 *CD ) / 1 8 . *R*RHOL*V*RH AT ) 

104700000 


C 



104800000 

003105 


64 

CONTINUE 

104900000 

003105 



IFI IBUGI7) . EQ. 0 ) RETURN 

105000000 

003107 



PRINT 66, I WII) ,1=1,14) , I GAMDOT I I ) , 1 = 1 , 15 ) , EM , TEM1 0,CDZ , TEM11 , CD, 

105100000 




1 VDOT 

105200000 

003143 


66 

FORMAT 1 4H W= / 2 I 6E20 . 8/ ) , 2 E2 C . 8 , 9H G AMD0T=/3 I 5E20 . 8/ ) / 29H EM, TEM 

105300000 




HO.CDZ.TEMl 1,CD,VD0T=/6E20.8/11H END BASIC/////) 

105400000 

003143 



RETURN 

105500000 


C 



105600000 

003144 



END 

105700000 


C 



105800000 


05 

CO 


05 

rf=- 


C00002 

000002 

C00002 

000002 

000002 

000002 

000002 

000002 


000002 

000004 

000012 

000013 

000016 

000020 

000020 

000026 

000026 

000027 

000030 

000032 

000033 


SUBROUTINE CHECK 

C *** SUBROUTINE CHECK, USED BY CALINT TO CONTROL SIZE OF COMPUTING 

C *** INTERVAL Cl 

C IN CHECK, CUVAR = UPCATED VALUES 

C VAR = VALUES AT BEGINNING OF INTERVAL 

C SPEC= COMPUTING INTERVAL USED ACROSS INTERVAL 

C Cl = COMPUTING INTERVAL FOR NEXT INTERVAL 

C IF ANSWERS NOT ACCEPTABLE, MODIFY SPEC, SET ELT=1 . AND RETURN 

C IF ANSWERS ACCEPTABLE, SET ELT=0 AND RETURN 

C 
C 

c 

COMMON VAR I 20) »DER(20),ELE1(20) ,ELE2(20),CI , SPEC , NDE , ELT ,NE RR 
COMMON RH3 L ,HL ,DN , E L , TW , C DC ,CDFM , FRO , CP VW , WASK , R, H A SK 

COMMON T 

C LABI IN MAIN, BASIC, FOFX CHECK 

COMMON/LABI /IBUG,CUVAR,A1,A2,A3,A4,A5,A6,TLIM1,TLIM2,TLIM3 
DIMENSION CUVAR ( 20 ) , IBUG120 ) 

C LAB2 IN MAIN, BASIC, CHECK, FOFX 

C0MM0N/LAB2 /ENTHAL 

C LAB3 IN MAIN , BASIC AND CHECK 

COMMON /LAB3/ TPRE V ,KEY , HEPS , KAY 
C 

DOUBLE PREC ISION VAR 
C 

C KEY SET =1 IN CHECK TO ALLOW TESTING OF ENTHALPY AFTER FIRST INTERVAL 

C KAY SET = 1 IN BASIC WHEN T.LT.O. 

C 

IFIKAY.EQ.l ) GO TO 10 

IF! ABSIT-TPREVJ .GT .500. ) GO TO 10 

IF(KEY.EQ.O) GO TO 20 

C CONTROL CHANGE IN ENTHAL 

IFtENTHAL. GT. 50000. ) GO TO 21 
IF( ENTHAL. LT. -500000. ) GO TO 21 
GO TO 20 

21 IF { AB S I (HPR EV-ENTHA L ) /ENTHAL ) .GT . HEPS ) GO TO 10 
20 CONTINUE 
TPREV=T 
ELT=0. 

HPREV=ENTHAL 

KEY=1 

RETURN 


105900000 
106000000 
106100000 
106200000 
106300000 
106400000 
106500000 
106600000 
106700000 
106800000 
106900000 
107000000 
107100000 
107200000 
107300000 
107400000 
107500000 
107600000 
107700000 
107800000 
107900000 
108000000 
108100000 
108200000 
108300000 
108400000 
108500000 
108600000 
108700000 
„ 108800000 
108900000 
109000000 
109100000 
109200000 
109300000 
109400000 
109500000 
109600000 
109700C00 
109800000 
109900000 
110000000 
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C 



110100000 

000034 


10 

ELT=1.0 

110200000 

000035 



IF { IBUG(IO) .EQ.O) GO TO 30 

110300000 

000037 



PRINT 11, T,TPREV,SPEC,ENTHAL,KEY,CI,HPREV 

110400000 

000060 


11 

FORMAT (4H T=E15.8, 8H TPREV=E15.8, 7H SPEC=E15.8, 9H ENTHAL= 

110500000 




1 E15.8, 6H K E Y= 1 1 / 32H REOUCE COMPUTING INTERVAL CI = E15.8/ 

110600000 




2 8H HPREV= El 5.8/) 

110700000 

000060 


30 

CONTINUE 

110800000 


C 


SPEC = COMPUTING INTERVAL JUST USED 

110900000 


C 


Cl = COMPUTING INTERVAL FOR NEXT STEP 

111000000 

000060 



SPEC=SPEC/2 . 

111100000 

000062 



CI=SPEC/2. 

111200000 

000063 


16 

CONTINUE 

111300000 

000063 



KEY = 0 

111400000 

000064 



KAY=0 

111500000 

000064 



IFICI .LT.l. 0E-20) PRINT 12 

111600000 

000073 


12 

FORMAT ( 50H SOMETHING MUST BE WRONG, C I . L T. 1 . 0E-20 , STOP 12/ 

111700000 




1 17H REEXAMINE INPUT) 

111800000 


C 


*** STOP 12 

111900000 

000073 



IFICI. LT.1.0E-20) STOP 12 

112000000 

000100 



RETURN 

112100000 

000101 



END 

112200000 


05 

C71 
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05 

a> 


000003 

000003 
000003 
000003 
C00003 
CO 0003 
000003 

000003 

000004 

000005 

C00006 

000007 


000010 

000011 

000014 

000026 

000040 


000043 

000072 

000073 

000113 


FUNCTION FOFX(TD) 

C FUNCTION FOFX(TD), USED BY ITR1 TO EVALUATE T 

C 

C LABI IN MAIN, BASIC, FOFX 

COMMON/L AB1 / IBUG,CUVAR,A1,A2,A3,A4,A5,A6,TLIM1,TLIM2,TLIM3 
C LAB2 IN MAIN, BASIC, CHECK, FOFX 

C0MM0N/LAB2 / ENTHAL 

DIMENSION All 15,4) ,A2( 15,4) , A3 ( 15,4! ,A4( 15,4) ,A5< 15,4 » ,A6< 15, 4) 
DIMENSION CUVAR120) , IBUGl 20 ) 

DIMENSION GAM( 15) 

EQUIVALENCE ( C UV AR< 4 ) , G AM ( 1 ) ) 

DATA REMD/8.31696E7/ 

C 

TB = TD 
T2=TD**2 
T3= T2*TD 
T4= T3*TD 
SUM=0. 

C TL I Ml = 1 COO . DEG.K 

C TL IM2=1 5000 » DEG K 

C TLIM3=1000000. 

DO 10 1=1,15 

IF(TD.LE.TLIMl) IT=1 

IF( TD.GT.TL IM1 .AND . TD. LE . TL IM2 ) IT=2 
IF! TD.GT.TL IM2 .AND. TD. LE.TLIM3) IT=3 
IFITD.GT.TLIM3) IT=4 
C 

C IH SUB I )/RT 

C 

HIORT = All I,IT)+A2(I,IT)/2. *TD + A3 ( I»IT)/3. *T2+A4< I, IT)/4.*T3 
1 +A5( I, IT ) / 5. *T4+A6 ( I, IT)/TD 
IF{ IBUGC8).. EQ.O) GO TO 10 
C 
C 

PRINT 20, I T,TD»HIORT ,GAMl I ) ,SUM 
10 SUM = SUM+ G AM ( I ) *H I ORT 

C R= 8 . 3 1 696E7 ERG/ MOLE-DEG 

C TC -ENTHAL/l 8 . 31 696E7*SUM ) 

C FOFX=TC 

C BECAUSE SUM SO SENSITIVE, TAKE IT OUT OF DENOMINATOR AND REDUCE 

C ITS EFFECT BY A FACTOR 9-8-67 

C EQUATION (12) 


112300000 

112400000 

112500000 

112600000 

112700000 

112800000 

112900000 

113000000 

113100000 

113200000 

113300000 

113400000 

113500000 

113600000 

113700000 

113800000 

113900000 

114000000 

114100000 

114200000 

114300000 

114400000 

114500000 

114600000 

114700000 

114800000 

114900000 

115000000 

115100000 

115200000 

115300000 

115400000 

115500000 

115600000 

115700000 

115800000 

115900000 

116000000 

116100000 

116200000 

116300000 

116400000 
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000122 


TC=TD+ . 1 *{ T D*REMD*SUM-ENTHAL ) 

116500000 

000127 


F0FX=TC 

116600000 


C 


116700000 

000130 


IF( IBUG(9) . EQ.O) RETURN 

116800000 


C 


116900000 


C 


117000000 

000131 


PRINT 20, TB,TC,SUM ,ENTHAL 

117100000 

000145 


20 FORMAT ( 5E1 8.8) 

117200000 


c 


117300000 

000145 


RETURN 

117400000 

000147 


END 

117500000 


oa 

-a 
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CT> 

oo 


000002 

000002 

000002 

000012 

000012 

000020 

000020 

000026 

000026 

000040 

000041 


SUBROUTINE ERROR 

C *** SUBROUTINE ERROR, USED BY CALINT IF INPUT IS BAD, BAD, BAD 

COMMON VAR< 20) ,DER< 20) ,ELE1(20) , EL E2 < 20 ) , C I , SPEC, NDE, EL T ,NERR 
C 

DOUBLE PRECISION VAR 
C 

PRINT 10, N ERR ,CI 

10 FORMAT ( 10H NERR= 15, 5H CI=E17.8/3 

IFINERR.EQ. 1) PRINT 11 

11 FORMAT ( 25H ***** BAD TNPUT *****/) 

IFINERR.EQ. 2) PRINT 12 

12 FORMAT I25H ***** ELE1=ELE2 *****/) 

C STOP 13 

IFINERR.EQ. 1. OR. NERR.EQ. 2) STOP 13 
RETURN 
END 


117600000 

117700000 

117800000 

117900000 

118000000 

118100000 

118200000 

118300000 

118400000 

118500000 

118600000 

118700000 

118800000 

118900000 

119000000 

119100000 
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SUBROUTINE CAL INK VAR, DER ,ELE1 , ELE2 ,C I , SPEC ,N, CUVAR ,ELT , 

119200000 




1CIMAX,NERR* PHMAX) 

119300000 


c 



119400000 


c 


CALINT, INTEGRATION ROUTINE 

119500000 


c 


THE NUMERICAL INTEGRATION ALGORITHM USED IS FOUND IN - A METHOD 

119600000 


c 


FOR THE NUMERICAL INTEGRATION OF COUPLED FIRST ORDER DIFFERENTIAL 

119700000 


c 


EQUATIONS WITH GREATLY DIFFERENT TIME CONSTANTS - BY CHARLES E. 

119800000 


c 


TREANOR, CONTRACT NO.NASR-119, CAL REPORT NO . AG-1 729-A-4, JANUARY 

119900000 


c 


1964. CORNELL AERONAUTICAL L ABORATORY , I NC CORNELL UNIVERSITY, 

120000000 


c 


BUFFALO, N.Y. 

120100000 


c 



120200000 

000017 



DOUBLE PRECISION VAR 

120300000 

000017 



DIMENSION VAR (42) ,DER(41) ,ELE1(41) , ELE2 ( 41 ) , CUV AR l 42) ,F1(41), 

120400000 




1F2I41 ) ,F3(41) ,CAPF1 (41) ,C APF 2 ( 41 ) , CAPF3 ( 4 1 ) , P ( 41 ) , PH ( 41 ) , 

120500000 




2DELTY (41) , Y 3 ( 4 1 ) , Y4 ( 41 ),F4(41 ) , Y2 ( 41 ) 

120600000 


c 



120700000 


c 


NERR= 1 Cl, N, OR, NT IS EQUAL TO ZERO 

120800000 


c 



120900000 


c 


NERR=2 ELE1 LESS THAN OR EQUAL TO ELE2 

121000000 


c 



121100000 


c 



121200000 


c 


TEST INPUT 

121300000 


c 



121400000 

000017 



FN=N 

121500000 

000020 



TEST=CI*FN 

121600000 

000021 



IF(TEST) 99 8,997, 998 

121700000 

000023 


997 

NERR= 1 

121800000 

000025 



GO TO 100 

121900000 

000025 


998 

IF ( ELE1-ELE 2) 999,999, 1000 

122000000 

000027 


999 

NERR=2 

122100000 

000031 



GO TO 100 

122200000 

000031 

1000 

IF! SPEC )5, 1,5 

122300000 


c 



122400000 


c 


SECTION FOR INITIALIZATION COMPUTATION OF DERIVATIVES 

122500000 


c 



122600000 

000032 


1 

SPEC=CI 

122700000 

000033 



ICGNT=1 

122800000 

000034 


2 

N1=N+1 

122900000 

000037 



DO 3 1=1, N1 

123000000 


c 



123200000 

000040 


3 

CUV AR ( I )=VAR( I ) 

123100000 

000047 



CALL BASIC 

123300000 


05 

CO 
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-3 

O 


C 

C RETURN WITH DERIVATIVES IN DER 

C 

000Q50 DO 4 1=1, N 

000055 4 F1(I)=0ER( I ) 

C 

000062 RETURN 

C 

C COMPUTE Y2,X2 

C 

000063 5 CUVAR(l)=VAR(l)+CI/2. 

000101 DO 6 1=1, N 

000102 11 = 1+1 

000103 Y2( I)=VAR( I D+CI/2. *F1(I) 

000127 IF(Y2 (I) )65 ,6,6 

000132 6 CUVARI I1) = Y2( I ) 

000141 GO TO 66 

000141 65 SPEC=C I 

000142 CI=CI/2. 

000143 GO TO 5 

C 

C CALL BASIC TO EVALUATE F2 

C 

000144 66 CALL BASIC 

C 

C RETURN 

C 

000145 DO 7 1=1, N 

000152 11=1+1 

000153 F2( II =DER( I ) 

C 

C COMPUTE Y3 

C 

000156 Y3 ( I ) =VAR( 11 ) +CI /2 . *F2 ( I ) 

000202 IF ( Y3 ( 1)165,7,7 

000205 7 CUVAR ( I1)=Y3( I ) 

C 

C CALL BASIC TO EVALUATE F3 

C 

000214 CALL BASIC 

C RETURN 

C 


000215 


DO 10 1=1, N 


123400000 
123500000 
123600000 
123700000 
123800000 
123900000 
124000000 
124100000 
124200000 
124300000 
124400000 
124500000 
124600000 
124700000 
124800000 
124900000 
125000000 
125100000 
125200000 
125300000 
125400000 > 

125500000 2 

125600000 M 

125700000 2 

125800000 O 

125900000 X 

126000000 
126100000 
126200000 
126300000 
126400000 
126500000 
126600000 
126700000 
126800000 
126900000 
127000000 
127100000 
127200000 
127300000 
127400000 
127500000 
127600000 



000222 



F3( I )=DER( I ) 

127700000 


C 



127800000 


c 


COMPUTE P* PH AND CAP F TERMS 

127900000 


c 



128000000 

000225 



IF(Y3( I) — Y 2 C I ) ) 9,8,9 

128100000 

000230 


8 

P(I)=0. 

128200000 

000232 



GO TO 91 

128300000 

000233 


q 

P(I > =— ( f F3 ( I)-F2(IJ )/(Y3( I)-Y2< HI) 

128400000 

000245 


91 

PH ( I )=P(I)*CI 

128500000 

000251 



IF ( PH ( I ) >83,83,103 

128600000 

000253 


83 

PH I I ) =0. 

128700000 

000255 



P(I)=0. 

128800000 

000257 



GO TO 84 

128900000 

000257 


103 

I F t ABS(Y3( I )-Y2( I) ) /( ( ABSIY3I I > ) + ABS(Y2(I > > ) /2 . >- . 5E-4 ) 83, 8 3 , 84 

129000000 

000272 


84 

IF ( PH { I )— . 1)85,85,95 

129100000 

000276 


85 

CAPF1 ( I ) =1 .-PHm/2. + (PH(I)**2)/6.-(PH( I)**3)/24. 

129200000 

000310 



CAPF2 ( I ) = . 5-PHI I)/6 . + { PHI I ) **2 ) /24.-I PHI I )**3)/120. 

129300000 

000323 



CAPF3 ( I )=l./6.-PH( I)/24.+(PH( I )**2)/120.-(PH( I )**3)/720. 

129400000 

000341 



GO TO 10 

129500000 

000342 


95 

CAP FI ( I ) = ( E XP t — PH ( I) ) — 1.)/<— PHI I) } 

129600000 

000352 



CAPF2 ( I ) = (CAPF1( I )-l.) /(-PHI I ) ) 

129700000 

000361 



CAPF3I I ) = (CAPF2( I )-.5> / C — PH ( I ) > 

129800000 

000371 


10 

CONTINUE 

129900000 


c 



130000000 


c 


IS PH BETWEEN ELE2 AND ELE1 

130100000 


c 



130200000 

000400 



IF ( I CONT-1 ) 101,101,102 

130300000 

000402 


102 

ICONT= ICONT-1 

130400000 

000404 



SPEC=CI 

130500000 

000405 



GO TO 17 

130600000 

000405 


101 

DO 11 1=1, N 

130700000 

000407 


104 

I F ( ABSIPHI I H-ELEH I ) ) 11, 11,13 

130800000 

000415 


11 

CONTINUE 

130900000 


c 



131000000 

000420 


12 

SPEC=CI 

131100000 

000421 



GO TO 15 

131200000 


c 



131300000 


c 


HALVE INTERVAL AND DOUBLE PH RANGE 

131400000 


c 



131500000 

000422 


13 

DO 96 1=1, N 

131600000 

000424 



ELE1( I )=ELE1( I >*2. 

131700000 

000426 



IFIELEl I I)-PHMAX)94, 94,955 

131800000 

000432 


94 

ELE2 ( I ) = ELE2( I )*2. 

131900000 
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-3 

ts3 


000435 GO TO 96 

000435 955 ELE1( I )=ELE1{ I )/2 . 

000440 96 CONTINUE 

000443 SPEC=C I 

000443 C I = C I / 2 . 

000444 ICONT=3 

000446 GO TO 5 

C 

C RETURN TO RECOMPUTE INTERVAL 

C 

c 

C RETURN TO RECOMPUTE INTERVAL 

C 

000446 15 DO 16 1=1, N 

000450 IF ( ABS ( PH ( I )! - EL E2 I I )>16 , 17 , 1 7 

000456 16 CONTINUE 

C 

C DOUBLE INTERVAL 

C 

000461 CI=2.*CI 

000462 I F I C I — C IMAX ) 17 ,17,165 

000465 165 C I = C I MAX 

000467 17 CONTINUE 

C 

C COMPUTE Y4,X4 

C 

000467 DO 18 1=1, N 

000471 11=1+1 

00 0472 CUV AR ( 1 1 ) = V AR ( I1) + SPEC*(F3( 1)4(2 

1 (CAPF1 IIJ-2 .*CAPF2t I ) ) +F2C I )*PH< 
000531 IFICUVARUl )) 175,18, 18 

000534 175 C 1 = SP EC 

000535 CI=CI/2. 

000537 GO TO 5 

000537 18 Y4( I ) =CUVAR (II) 

C 

000546 CUV AR ( 1 )=VAR( 1 J+SPEC 

C 

C CALL BASIC TO EVALUATE F4 

C 

000553 CALL BASIC 

C RETURN 

C 


*CAPF2<I))+F1( I ) * 
)*CAPF2( I) I 


132000000 
132100000 
132200000 
132300000 
132400000 
132500000 
132600000 
132700000 
132800000 
132900000 
133000000 
133100000 
133200000 
133300000 
133400000 
133500000 
133600000 
133700000 
133800000 
133900000 
134000000 
134100000 
134200000 
134300000 
134400000 
134500000 
134600000 
134700000 
134800000 
134900000 
135000000 
135100000 
135200000 
1 35300000 
135400000 
135500000 
135600000 
135700000 
135800000 
135900000 
136000000 
136100000 
136200000 
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000554 


DO 20 1=1, N 

136300000 



000561 


11=1+1 

136400000 



C00562 


F4( I ) =DER( I ) 

136500000 




C 


136600000 




C 

COMPUTE DELTA Y 

136700000 




C 


136800000 



000565 


DELTY ( I )=SPEC*(F1< I )*CAPF1( I )+(-3.*{Fl( I) +P ( I )*VAR( 11 ) ) 

136900000 





1+2.*(F2 ( I)+P( I )*Y2( I) ) +2 . * ( F3 ( I ) + P ( I ) *Y3 ( I ) )-F4 ( I ) 

137000000 





2-P(I)*Y4( I) )*CAPF2<I)+4.*((F1(I)+P( I >*VAR(I1I)-(F2( I) 

137100000 





3 + P( I ) *Y2 ( I ) )- ( F3 ( I } +P ( I ) *Y3 ( I )) + {F4( I )+P( I ) *Y4 ( I) ) )* 

137200000 





4C APF3 ( I ) ) 

137300000 




C 


137400000 




C 

COMPUTE Y+DELTA Y 

137500000 




c 


137600000 



000742 


20 CUVAR ( I1) = VAR( I1)+0ELTY(I) 

137700000 




c 


137800000 




c 

CALL CHECK FOR DECISION TO ACCEPT OR RECOMPUTE 

137900000 




c 

INTERVAL 

138000000 



000757 


CALL CHECK 

138100000 

t> 


000760 


IF ( ELT )21 , 21 , 23 

138200000 




c 


138300000 

►d 



c 

UPDATE Y VALUES 

138400000 




c 


138500000 



000766 


21 N 1=N+ 1 

138600000 

d 


000771 


DO 22 1 = 2, N1 

138700000 



000772 


11=1-1 

138800000 



000773 


22 VAR(I)=VAR( D+DELTYU1) 

138900000 



001006 


VAR ( 1 ) =VAR ( 1 ) + SPEC 

139000000 




c 


139100000 




c 

RETURN TO COMPUTE DERIVATIVES AT Y+DELTA Y 

139200000 




c 


139300000 



001014 


GO TO 2 

139400000 




c 


139500000 




c 

RETURN TO RECOMPUTE INTERVAL 

139600000 




c 


139700000 



001014 


23 GO TO 5 

139800000 




c 


139900000 



001015 


100 CALL ERROR 

140000000 



001016 


RETURN 

140100000 



001017 


END 

140200000 








CO 
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TABLE II.- CONSTANTS FOR USE IN REACTION 
RATE FORMULAS 


j Reaction 

1 h 2 o + M xOH + H + M 

2 0 2 + M^0+0 + M 

3 N 2 + M= N+ N + M 

4 H 2 + M:H + H + M 

5 NO + M= N + 0 + M 

6 OH + MsO + H + M 

7 H + 0 2 X OH + O 

8 O + H 2 — OH + H 

9 H + NO X OH + N 

10 0+N0:0 2 +N 

11 O +N 2 r NO +N 

12 O + H 2 0 = OH + OH 

13 H + H 2 0 sOH + H 2 

14 NO + H 2 x H 2 0 + N 

15 N + 0;N0 + +e- 

16 N + N ~ N 2 + + e“ 

17 N 2 + e" Ng” 1 " + e" + e“ 

18 NO + e" r NO + + e“ + e‘ 

19 N + e" X N + + e“ + e” 

20 O + e" X 0 + + e“ + e" 

21 H + e~ r H + + e“ + e" 


/J-j, g/mol 

10.957 

14.933 

14.000 
1.867 

14.483 

10.578 

.970 

1.778 

.968 

10.435 

10.182 

8.471 

.947 

1.875 

7.467 

7.000 

5.486 X lO" 4 
5.486 X 10“ 4 
5.486 X 10~ 4 
5.486 X 10 -4 
5.486 X 10 -4 


Dj oj, cm^ 

217.10 10.00 X 10“ 16 

217.34 10.00 

414.50 10.00 

190.34 10.00 

275.74 10.00 

186.09 10.00 

31.44 1.00 

4.25 1.00 

89.65 1.00 

58.21 1.00 

138.93 2.69 

31.02 1.00 

26.77 1.00 

62.88 1.00 

117.69 .018 

247.27 .066 

661.94 .879 

393.00 .879 

617.75 .879 

578.24 .879 

577.81 .879 
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TABLE II.- CONSTANTS FOR USE IN REACTION 


RATE FORMULAS - Concluded 


j 

Reaction 

Mj, g/mol 

D i 


cm^ 

22 

N 2 + NO + - N 2 + + NO 

14.483 

268.94 

20.00 X 10“ 16 

23 

N 2 + N + = N 2 + + N 

9.333 

44.19 



24 

N 2 + 0 + 5 : N 2 + + O 

10.182 

83.70 



25 

N 2 + H + X N 2 + + H 

.966 

84.12 



26 

N + NO+ X N + + NO 

9.545 

224.75 



27 

O + NO+ tO + +NO 

10.435 

185.24 



28 

H + NO + - H + + NO 

.968 

184.82 



29 

N + 0 + X N + + O 

7.467 
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Figure 4.- Velocity along body streamline at an altitude 
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47.55 km. (Except as noted, results obtained by using curve A of fig. 1(b).) 





Figure 5.- Fraction of water evaporated along body streamline at an altitude of 47.55 km. (Except as noted, results obtained by using curve A of fig. 1(b).) 
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